Last summer, while relaxing on the beaches of Berck, a French town known for treating tuberculosis in kids by exposing them to the fresh maritime air (back in 19th century, they have antibiotics these days), I found myself daydreaming about building my own programming language.

Spoiler alert: I don’t know how to build programming languages, but I have developed extremely strong opinions over the years about the features a modern data science language should have. So could I use them fancy LLMs to build one?

Also, let’s get one question answered straight away: why create a new language instead of contributing to existing ones? I certainly do contribute, I maintain several R packages like {rix}, {rixpress}, and {chronicler}, and even have two Python packages (cronista and ryxpress), but I wanted a clean slate to build a system centered around a few non-negotiable principles and features I’ve implemented over the years in R:

• Reproducibility-First: A language where reproducibility isn’t a bolt-on afterthought managed by external tools, but the very foundation of the runtime.
• Aggressive Re-use: Instead of reinventing the wheel, this language would stand on the shoulders of giants. It’d use Nix for package management and environment isolation, and Apache Arrow as its high-performance backbone for data frames. R, Python, Julia and other languages would provide the algorithms and models.
• First-Class Pipelines: Scripts shouldn’t be a sequence of side-effects. In this language, pipelines would be mandatory and first-class citizens.
• Fail Early and Loudly: No silent type conversions or hidden NAs. If something is wrong, the language breaks immediately so you can fix it.
• Errors as Objects: Inspired by functional programming, errors are first-class values that can be inspected and handled gracefully.
• Two Pipes: I want two pipes, one for linear transformations, |>, and a maybe-pipe, ?|> for error recovery. Unlike the standard pipe, ?|> always forwards its value, including Errors, to the next function, allowing you to write handlers that inspect and potentially recover from them. Since Errors are just values, this composes naturally with the rest of the language.
• Polyglot by Design: Rather than re-implementing every statistical algorithm, this language would be designed to orchestrate and bridge R, Python, and Julia seamlessly.

Also, we’re in a post LLM world, and like them or not, they’re here to stay. They’re pretty useful to write boilerplate code and so any new language would be dead on arrival if it didn’t play nicely with LLMs. So such a new language would need to be written for LLMs primarily, because I don’t expect anyone to learn any new language. This is where the declarative nature of Nix is a huge advantage. Because environments are precisely described, it is much easier for LLMs to focus on generating code and not have to fight with environment setup. This is also the reason I took another radical decision: since Nix would be mandatory for setting up the environment, why bother building OS-specific binaries? I’d just build a Nix package for this language and let Nix handle the rest.

This architecture results in a DSL for orchestration, making it trivial to transfer data objects between different ecosystems without the usual FFI (Foreign Function Interface) friction.

With these ideas in mind, I started prompting Gemini to brainstorm and started by generating specification files. Very broad first, but as days went by, more and more focused. The way I went about it (and still go) is that I first brainstorm an idea with an LLM, then I ask it to generate a specification file, then I refine it, ask it to generate a new specification file, and so on. Once I’m happy with the spec, I ask an LLM to generate a minimal implementation of the spec. Usually writing the spec and a first implementation is a task shared between Claude and Gemini (through Antigravity). Then I open a pull request and ask GitHub Copilot to review it (usually with GPT-5.x). I repeat this process until I’m happy with the implementation. I always ask for documentation and unit tests (and golden tests when relevant, more on this later).

I started to really believe that I had something interesting, so I gave it a shot, and called it T. I had long joked that the natural successor to R should be called T (because R is the successor to S… and no, I’m not going to call it Q because that sounds like the word for ass in French).

Something else that made me confident I could succeed, besides my own hubris, was that I am pretty familiar with unit testing, test-driven development, trunk-based development and Nix. When you combine all these elements, it makes developing with LLMs quite safe.

So I just started prompting. And now I’m quite happy to announce that there is a beta version of T that you can use today!

By leveraging Nix as a build engine, T can treat complex data science workflows as buildable derivations. A typical T pipeline looks like this:

p = pipeline {
  -- 1. Python node: read data with pandas
  mtcars_pl = pyn(
    command = <{
import pandas as pd
pd.read_csv("data/mtcars.csv", sep="|")
    }>,
    include = ["data/mtcars.csv"],
    serializer = ^csv
  )

  -- 2. Python node: filter and serialize as CSV
  mtcars_pl_am = pyn(
    command = <{
mtcars_pl[mtcars_pl['am'] == 1]
    }>,
    deserializer = ^csv,
    serializer = ^csv
  )

  -- 3. R node: read CSV and take head using functions.R
  mtcars_head = rn(
    command = <{
my_head(mtcars_pl_am)
    }>,
    functions = ["src/functions.R"],
    deserializer = ^csv,
    serializer = ^csv
  )

  -- 4. R node: select column with dplyr
  mtcars_mpg = rn(
    command = <{
library(dplyr)
mtcars_head %>% select(mpg)
    }>,
    deserializer = ^csv,
    serializer = ^csv
  )

  -- Render Quarto report
  report = node(script = "src/report.qmd", runtime = Quarto)
}

-- Materialize the pipeline
populate_pipeline(p, build = true)
pipeline_copy() -- Copy the outputs from the Nix store to your working directory

As you can see, each node has a command argument where you can write literal R or Python code. It is also possible to provide the path to a script instead. If packages need to be loaded for the code to work, you can just write the calls to load the required packages in the command argument as well.

While T is heavily inspired by the {targets} package in R, it takes the concept a step further by making pipelines first-class objects within the language itself. This means you can:

• Compose Pipelines: You can define small, modular pipelines and then merge them into larger ones using standard operators.
• Static Analysis: Because the DAG (Directed Acyclic Graph) is defined within the language, T can validate your entire workflow (checking for circular dependencies or missing data) before a single line of code even runs.
• Heterogeneous Execution: A single pipeline can effortlessly mix R, Python, and native T code. Data is passed between these nodes using built-in serializers like ^csv, ^arrow, or even specialized formats like ^pmml for traditional models and ^onnx for deep learning architectures. It is also possible to define your own serializers.
• Immutable State: Each node output is managed by Nix, meaning if you haven’t changed the code or the data for a specific node, T (via Nix) will simply pull the cached result from previous runs.

But don’t let the “orchestrator” label fool you; T is also a capable language in its own right. It features a selection of built-in packages inspired by the tidyverse for data manipulation. Thanks to its Arrow backend, it is surprisingly fast. I even maintain a CI benchmark running on NYC Taxi data to ensure performance remains competitive.

I made sure that T is pretty easy to use with LLMs by providing a file called summary.md in the root of the GitHub repository. This file is meant to be used by LLMs to quickly learn the language’s syntax and generate code accordingly. You could also provide the whole help documentation to the LLM (found in the repository under help/docs.json), but I found that a summary is usually enough. There is also another experimental feature I’m thinking about, called intent blocks. These blocks would essentially be first-class structured comments that would be used to anchor LLM’s behaviour and make it more deterministic. These blocks would be parsed by T and used to generate code accordingly. I have some ideas how these could look like, something like this:

intent {
  description: "Customer churn prediction",
  assumptions: ["Age > 18", "NA imputed with multiple imputation"],
  requires: ["dataset.csv"]
}

# 1. Initialize a new project
nix run github:b-rodrigues/tlang -- init --project my_t_project
cd my_t_project

[project]
name = "r_py_xgboost_t"
description = "A T data analysis project"

[dependencies]
# T packages this project depends on
# Format: package = { git = "repository-url", tag = "version" }
# Example:
# stats = { git = "https://github.com/t-lang/stats", tag = "v0.5.0" }

[r-dependencies]
packages = ["dplyr", "yardstick"]

[py-dependencies]
version = "python313"
packages = ["numpy", "pandas", "scikit-learn", "xgboost"]

[additional-tools]
packages = ["quarto"]

[t]
# Minimum T language version required
min_version = "0.51.2"

\(\)
packages = ["amsmath", "geometry", "hyperref", "biblatex"]

nix run github:b-rodrigues/tlang -- init --package my_package
cd my_package

[package]
name = "my_package"
version = "0.1.0"
description = "A brief description of what my_package does"
authors = ["brodriguesco"]
license = "EUPL-1.2"
homepage = ""
repository = ""

[dependencies]
# T packages this package depends on
# Format: package = { git = "repository-url", tag = "version" }

[t]
# Minimum T language version required
min_version = "0.5.0"

nix develop

Agentic AI tools like Claude Code can write and run code, fix its own errors, and produce a formatted report with figures. I wanted to know whether that translates into reliable ecological modelling, so we ran a test: three fisheries tasks, four AI models, ten independent runs each, scored against a rubric. The results are published in Fish and Fisheries.

We found agents can be genuinely useful, but only if you know how to use them well and only if you know enough about the analysis to catch what they miss.

A little bit less than two years ago, building on work by Jim Hester and Kevin Ushey, Davis Vaughan completed a very impactful JavaScript file for the R community: an R grammar for the Tree-sitter parsing generator. He even got a round of applause for it during a talk at the useR! 2024 conference! So, did he get cheered for… grammatical rules in a JavaScript file?

No, the audience was excited about the improved developer experience for R that this file unlocked. R tooling around Tree-sitter is how you get

• reformatting through Air and linting through Jarl;
• auto-completion or help on hover in the Positron IDE;
• better search for R on GitHub;
• and more!

In this post, we’ll explain what Tree-sitter is, and how tools built on Tree-sitter can benefit your R development workflow.

Code parsing: what is Tree-sitter?

Tree-sitter is a code parsing generator written in C, with bindings existing in several languages including Rust (and R!).

Let’s rewind a little bit. What does it mean to parse code?

Basically, given a string of code like

a <- mean(x, na.rm = TRUE)

How do you know that mean is a function name, na.rm an argument name, TRUE a logical? You have to parse that code into what’s called a parse tree. You do that in your head when reading R code.

R itself can obviously parse R code, thanks to its grammar. See for instance the commit that introduced R’s native pipe, which necessitated extending R’s syntax thus modifying its grammar.

You can use parse() and getParseData() to parse R code.

parse(
 text = "a <- mean(x, na.rm = TRUE)",
 keep.source = TRUE
) |>
 getParseData()
#> line1 col1 line2 col2 id parent token terminal text
#> 23 1 1 1 26 23 0 expr FALSE 
#> 1 1 1 1 1 1 3 SYMBOL TRUE a
#> 3 1 1 1 1 3 23 expr FALSE 
#> 2 1 3 1 4 2 23 LEFT_ASSIGN TRUE <-
#> 21 1 6 1 26 21 23 expr FALSE 
#> 4 1 6 1 9 4 6 SYMBOL_FUNCTION_CALL TRUE mean
#> 6 1 6 1 9 6 21 expr FALSE 
#> 5 1 10 1 10 5 21 '(' TRUE (
#> 7 1 11 1 11 7 9 SYMBOL TRUE x
#> 9 1 11 1 11 9 21 expr FALSE 
#> 8 1 12 1 12 8 21 ',' TRUE ,
#> 13 1 14 1 18 13 21 SYMBOL_SUB TRUE na.rm
#> 14 1 20 1 20 14 21 EQ_SUB TRUE =
#> 15 1 22 1 25 15 16 NUM_CONST TRUE TRUE
#> 16 1 22 1 25 16 21 expr FALSE 
#> 17 1 26 1 26 17 21 ')' TRUE )

Or you could transform that same data into XML using Gábor Csárdi’s {xmlparsedata}:

parse(
 text = "a <- mean(x, na.rm = TRUE)",
 keep.source = TRUE
) |>
 xmlparsedata::xml_parse_data(pretty = TRUE) |>
 xml2::read_xml() |>
 as.character() |>
 cat()
#> <?xml version="1.0" encoding="UTF-8" standalone="yes"?>
#> <exprlist>
#> <expr line1="1" col1="1" line2="1" col2="26" start="28" end="53">
#> <expr line1="1" col1="1" line2="1" col2="1" start="28" end="28">
#> <SYMBOL line1="1" col1="1" line2="1" col2="1" start="28" end="28">a</SYMBOL>
#> </expr>
#> <LEFT_ASSIGN line1="1" col1="3" line2="1" col2="4" start="30" end="31"><-</LEFT_ASSIGN>
#> <expr line1="1" col1="6" line2="1" col2="26" start="33" end="53">
#> <expr line1="1" col1="6" line2="1" col2="9" start="33" end="36">
#> <SYMBOL_FUNCTION_CALL line1="1" col1="6" line2="1" col2="9" start="33" end="36">mean</SYMBOL_FUNCTION_CALL>
#> </expr>
#> <OP-LEFT-PAREN line1="1" col1="10" line2="1" col2="10" start="37" end="37">(</OP-LEFT-PAREN>
#> <expr line1="1" col1="11" line2="1" col2="11" start="38" end="38">
#> <SYMBOL line1="1" col1="11" line2="1" col2="11" start="38" end="38">x</SYMBOL>
#> </expr>
#> <OP-COMMA line1="1" col1="12" line2="1" col2="12" start="39" end="39">,</OP-COMMA>
#> <SYMBOL_SUB line1="1" col1="14" line2="1" col2="18" start="41" end="45">na.rm</SYMBOL_SUB>
#> <EQ_SUB line1="1" col1="20" line2="1" col2="20" start="47" end="47">=</EQ_SUB>
#> <expr line1="1" col1="22" line2="1" col2="25" start="49" end="52">
#> <NUM_CONST line1="1" col1="22" line2="1" col2="25" start="49" end="52">TRUE</NUM_CONST>
#> </expr>
#> <OP-RIGHT-PAREN line1="1" col1="26" line2="1" col2="26" start="53" end="53">)</OP-RIGHT-PAREN>
#> </expr>
#> </expr>
#> </exprlist>

In both cases, you recognize words such as LEFT_ASSIGN or SYMBOL_FUNCTION_CALL. Parsing is an essential step before the code is actually executed, but parsed code can also be used for other purposes, such as analyzing code without brittle regular expressions (does it call a particular function?), navigating code (going from a function call to the definition of that function), or modifying code (replacing all occurrences of a function with another one).

Now, Tree-sitter performs this same code parsing but faster especially thanks to its support of incremental parsing – which is key to updating the syntax tree as you are typing in your editor for instance! Tree-sitter is agnostic in that it can parse any code as long as there is a grammar for it (think, Rosetta Stone plugins). It’s been used for many languages which means many tools have been built around it.

To have Tree-sitter “learn” a new language you need to give it a file containing the definition of the syntax of that language, what’s called a grammar. This is where the aforementioned JavaScript file by Davis Vaughan and collaborators comes into play! The treesitter-r repo, which provides a translation of the R grammar in the format expected by Tree-sitter, is the base of all tools presented in this post which use R code as their input.

Here’s how to use the {treesitter} R package for the same code as earlier. The {treesitter} R package allows us to use Tree-sitter from R. To parse R code with it, we need the language() function from {treesitter.r}¹.

library(treesitter)
#> 
#> Attaching package: 'treesitter'
#> The following object is masked from 'package:base':
#> 
#> range
language <- treesitter.r::language()
parser <- parser(language)
text <- "a <- mean(x, na.rm = TRUE)"
parser_parse(parser, text)
#> <tree_sitter_tree>
#> 
#> ── Text ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#> a <- mean(x, na.rm = TRUE)
#> 
#> ── S-Expression ────────────────────────────────────────────────────────────────────────────────────────────────────────────────
#> (program [(0, 0), (0, 26)]
#> (binary_operator [(0, 0), (0, 26)]
#> lhs: (identifier [(0, 0), (0, 1)])
#> operator: "<-" [(0, 2), (0, 4)]
#> rhs: (call [(0, 5), (0, 26)]
#> function: (identifier [(0, 5), (0, 9)])
#> arguments: (arguments [(0, 9), (0, 26)]
#> open: "(" [(0, 9), (0, 10)]
#> argument: (argument [(0, 10), (0, 11)]
#> value: (identifier [(0, 10), (0, 11)])
#> )
#> (comma [(0, 11), (0, 12)])
#> argument: (argument [(0, 13), (0, 25)]
#> name: (identifier [(0, 13), (0, 18)])
#> "=" [(0, 19), (0, 20)]
#> value: (true [(0, 21), (0, 25)])
#> )
#> close: ")" [(0, 25), (0, 26)]
#> )
#> )
#> )
#> )

Tree-sitter is the workhorse of many tools, that are mentioned in the diagram below. All of them are dependent on tree-sitter and the R grammar provided to it. Some of them are command-line interfaces (CLIs), while others are R packages.

Browsing code interactively: Positron IDE, GitHub

The real reason why the audience applauded Davis Vaughan is that he explained how the R grammar for Tree-sitter had been deployed to GitHub so that we get almost as good experience browsing R code on GitHub as browsing, say, JS code. If we search for a function name in a repository for instance, its definition will be indicated in the search results. See Davis’ slides (also available in PDF), or refer to the video below showing how typing vetiver_model in the search bar from the R vetiver repo makes the function definition the first result, on which one can click to land into the definition.

Also very useful is the use of Tree-sitter by Ark, the R kernel used in the Positron IDE. Ark is how you get autocompletion and help on hover in Positron. The video below shows how you can extend the selection to further steps of a pipeline in Positron.

This use case of Tree-sitter is also featured in Davis’ slides. See also Lionel Henry’s and Davis Vaughan’s talk about Ark at posit conf 2024, especially the part about code assistance.

Other development environments such as Emacs also have support for Tree-sitter.

Searching/browsing code

You can parse and search R code using the {treesitter} R package and treesitter query syntax. The {treesitter} R package is a dependency of the {gander} package by Simon Couch, that is meant to be used for a better experience with LLMs when writing R code. Another use case of the {treesitter} R package is the {igraph.r2cdocs} extension to {roxygen2} for the {igraph} package, that parses all of igraph R code to then be able to identify, for each exported function, whether it (in)directly calls a function whose name ends with _impl, indicating a wrapper to a C igraph function whose docs can be then be linked from the manual of the R function.

The {pkgdepends} package calls Tree-sitter (C) to detect dependencies in files. Below we run it on the source of the saperlipopette R package.

pkgdepends::scan_deps(
 "../../../../../CHAMPIONS/saperlipopette",
 "../../../../../CHAMPIONS"
)
#> 
#> Dependencies:
#> + brio  @ R/blame.R, R/check-editor.R, R/clean-dir.R, R/committed-to-main.R, R/committed-to-wrong-branch.R, R/conflict…
#> + cli  @ inst/exo_bisect-Rprofile.en.R, inst/exo_bisect-Rprofile.es.R, inst/exo_bisect-Rprofile.fr.R, inst/exo_blame-…
#> + devtools  @ saperlipopette.Rproj
#> + fs  @ R/blame.R, R/check-editor.R, R/clean-dir.R, R/committed-to-main.R, R/committed-to-wrong-branch.R, R/conflict…
#> + gert  @ inst/exo_check_editor-Rprofile.en.R, inst/exo_check_editor-Rprofile.es.R, inst/exo_check_editor-Rprofile.fr.…
#> + knitr  @ README.Rmd
#> + parsedate  @ R/utils-git.R
#> + purrr  @ R/create-all.R, R/debug.R, R/log-deleted-file.R, R/log-deleted-line.R, R/revparse.R, R/roxygen2.R, R/worktre…
#> + rlang  @ R/create-all.R, R/roxygen2.R, R/utils-fs.R, R/utils-usethis.R, R/zzz.R
#> + rmarkdown  @ README.Rmd, vignettes/saperlipopette.qmd
#> + roxygen2  @ R/roxygen2.R, saperlipopette.Rproj
#> + saperlipopette @ README.Rmd, vignettes/saperlipopette.qmd
#> + tibble  @ R/roxygen2.R
#> + usethis  @ R/blame.R, R/check-editor.R, R/clean-dir.R, R/committed-to-main.R, R/committed-to-wrong-branch.R, R/conflict…
#> + vctrs  @ R/roxygen2.R
#> + withr  @ R/blame.R, R/check-editor.R, R/clean-dir.R, R/committed-to-main.R, R/committed-to-wrong-branch.R, R/conflict…
#> 
#> Test dependencies:
#> + fs  @ tests/testthat/test-blame.R, tests/testthat/test-check-editor.R, tests/testthat/test-clean-dir.R, tests/test…
#> + gert  @ tests/testthat/test-blame.R, tests/testthat/test-clean-dir.R, tests/testthat/test-committed-to-main.R, tests…
#> + rlang  @ tests/testthat/test-blame.R, tests/testthat/test-check-editor.R, tests/testthat/test-clean-dir.R, tests/test…
#> + saperlipopette @ tests/testthat.R
#> + testthat  @ tests/testthat.R
#> + withr  @ tests/testthat/test-blame.R, tests/testthat/test-check-editor.R, tests/testthat/test-clean-dir.R, tests/test…

ast-grep is a useful tool built on Tree-sitter for searching and re-writing code, with a clearer query syntax than Tree-sitter’s. Its name is reminiscent of grep, but with ast-grep we do not need to write brittle regular expressions . {astgrepr} by Etienne Bacher is an R wrapper to the Rust bindings of ast-grep, and is used in Etienne’s {flir} package for refactoring code.

The ast-grep command-line interface (CLI) itself is featured in a useful blog post by Emil Hvitfeldt where he explains how to document the usage of ast-grep for Claude.

Formatting and linting: Air, Jarl

Speaking of CLIs…

Air, by Davis Vaughan and Lionel Henry, is a CLI built on Tree-sitter, in Rust. It reformats code blazingly fast.

Jarl, by Etienne Bacher, is a CLI built on Air, therefore also on Tree-sitter, in Rust. It lints and fixes code, also blazingly fast. It can even detect unreachable code, unused functions and duplicated function definitions.

In both of these examples, the creation of CLIs wrapping Rust bindings was more efficient than the creation of R packages wrapping the {treesitter} R package, for several reasons:

• Rust CLIs can edit code very fast²;
• CLIs are integrated in extensions for popular IDEs (for instance Positron);
• a CLI is easier to install on CI than an R package that needs, well, an R installation.

More tools

A brief mention of some other interesting tools we’ve explored a bit less.

Configuring: {ts} for parsing JSON and TOML (not R!)

The {ts} package by Gábor Csárdi is the backbone of two R packages used for editing and manipulating:

• TOML {tstoml};
• JSON {tsjson}.

Compared to existing parsers in R for those formats, these two packages preserve comments.

Testing code: {muttest}

Mutation testing is a kind of testing where you, say, randomly swap + with - in your code (you mutate it) and you run your tests to see whether they catch the mutant. The {muttest} package by Jakub Sobolewski is an R package for mutation testing, that depends on the {treesitter} R package.

Diffing code: difftastic

The difftastic CLI by Wilfred Hughes is “a structural diff tool that understands syntax”. This means that difftastic doesn’t only compare line or “words” but actual syntax by looking at lines around the lines that changed (by default, 3). Even better, it understands R out of the box. See this blog post with examples of R code diffing.

Conclusion: more to come?

In this post, we’ve presented an overview of Tree-sitter based tooling for R or in R.

Note that this ecosystem of tools is very actively developed, so some tools might come and go. However, the idea that plugging the R grammar into a general parsing generator brings cool features to us R developers, will remain true. Maybe you will contribute to this ecosystem, either through an existing tool or by creating a new one?

IMPORTANT: The website https://www.techtonique.net is down until further notice.

https://www.techtonique.net contained an language-agnostic API for machine learning tasks (classification, regression, survival analysis, forecasting etc.).

As a result, do not buy the Gumroad tutorial then.

You can still use the packages https://github.com/Techtonique locally.

PS: It’s not an April’s fool joke.

Disclaimer: This blog contains opinions that are of the authors alone and do not necessarily reflect the strategy of their respective organizations.

install.packages("dataviewR")

SEX == "F" & AGE > 65 & TRT01P == "Xanomeline High Dose"

library(dataviewR)
library(pharmaverseadam)

dataviewer(adsl, adlb)

@online{pujari2026,
  author = {Pujari, Siddhesh and Kumar N, Madhan and S, Gomathi and
    Haight, Mackenzie},
  title = {Meet {dataviewR:} {The} {View()} {You} {Always} {Wanted}},
  date = {2026-03-31},
  url = {https://pharmaverse.github.io/blog/posts/2026-03-29-meet-dataviewr-the/meet-dataviewr-the-view-you-always-wanted.html},
  langid = {en}
}

The Emergency departments leading the transformation of Alzheimer’s and dementia care (ED-LEAD) study, which I have written about in the past, is approaching the end of its third year. This multifactorial design evaluates three independent, yet potentially synergistic, interventions aimed at improving care for persons living with dementia (PLWD) and their caregivers.

To estimate intervention effects, we are using what I’ve called the HEx-factor model, a Bayesian hierarchical exchangeable factorial model. The original plan was to conduct all analyses using Stan. However, we’ve run into a bit of a snafu. I’ve been working through the problem, and thought I’d share here.

The challenge turns out to be a computational one. Because the ED-LEAD analyses must be conducted on National Institute on Aging (NIA) Data LINKAGE servers, we are working in a somewhat restricted software environment, at least with respect to Bayesian data analysis. In particular, we have not been able to install or run Stan, which was our analytic engine of choice. This forced us to consider alternatives, and we turned to JAGS, which is available in the Linkage environment and certainly is well-suited for Bayesian hierarchical modeling.

At first glance, this might seem like a straightforward substitution. Both Stan and JAGS allow us to specify the same likelihood and priors. However, I quickly noticed that the models were not performing as well in JAGS as they had in Stan. It turns out that the samplers used in JAGS are more sensitive to posterior dependence than the Hamiltonian Monte Carlo (HMC) methods implemented in Stan.

I set out to understand and fix the problem, and found that a simple reparameterization—re-coding the binary treatment indicators—made a substantial difference. With this change, the JAGS sampler was able to explore the posterior distribution much more efficiently, yielding results comparable to those obtained with Stan.

To understand why this happens, I ran a series of simple simulations comparing the original and reparameterized versions of a basic two-way factorial model. That is what I present here.

library(simstudy)
library(data.table)
library(ggplot2)
library(rjags)
library(coda)
library(posterior)
library(broom)
library(gt)

RNGkind("Mersenne-Twister", "Inversion", "Rejection")
set.seed(824)

s_gen <- function(n = 2000,
                    alpha = -0.8,
                    beta_a = 0.5,
                    beta_b = 0.9,
                    beta_ab = -0.3) {
  
  def <- 
    defData(varname = "A", formula = 0.5, dist = "binary") |>
    defData(varname = "B", formula = 0.5, dist = "binary") |>
    defData(varname = "AB", formula = "A*B", dist = "nonrandom") |>
    defData(varname = "A_c", formula = "A - 0.5", dist = "nonrandom") |>
    defData(varname = "B_c", formula = "B - 0.5", dist = "nonrandom") |>
    defData(varname = "AB_c", formula = "A_c * B_c", dist = "nonrandom") |>
    defData(
      varname = "Y", 
      formula = "..alpha + ..beta_a * A + ..beta_b * B + ..beta_ab * AB",
      dist = "binary", link = "logit"
    )
    
  genData(n, def)
  
}

dd <- s_gen()

fit_01 <- glm(Y ~ A * B, data = dd, family = binomial)
fit_c  <- glm(Y ~ A_c * B_c, data = dd, family = binomial)

tidy(fit_01)
## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -1.03     0.0995    -10.4  3.90e-25
## 2 A              0.835    0.135       6.18 6.34e-10
## 3 B              1.14     0.134       8.46 2.73e-17
## 4 A:B           -0.670    0.186      -3.61 3.11e- 4
tidy(fit_c)
## # A tibble: 4 × 5
##   term        estimate std.error statistic  p.value
##   <chr>          <dbl>     <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -0.212    0.0464     -4.57 4.91e- 6
## 2 A_c            0.501    0.0929      5.39 7.05e- 8
## 3 B_c            0.802    0.0929      8.63 6.07e-18
## 4 A_c:B_c       -0.670    0.186      -3.61 3.11e- 4

model_01 <- "
model {
  for (i in 1:N) {
    Y[i] ~ dbern(p[i])
    logit(p[i]) <- alpha + beta_a * A[i] + beta_b * B[i] + beta_ab * AB[i]
  }
  
  alpha   ~ dnorm(0, 0.25)
  beta_a  ~ dnorm(0, 0.25)
  beta_b  ~ dnorm(0, 0.25)
  beta_ab ~ dnorm(0, 25)
}
"

model_c <- "
model {
  for (i in 1:N) {
    Y[i] ~ dbern(p[i])
    logit(p[i]) <- alpha + gamma_a * A_c[i] + gamma_b * B_c[i] + gamma_ab * AB_c[i]
  }
  
  alpha   ~ dnorm(0, 0.25)
  gamma_a  ~ dnorm(0, 0.25)
  gamma_b  ~ dnorm(0, 0.25)
  gamma_ab ~ dnorm(0, 25)
}
"

fit_jags <- function(dat, model_string, centered = FALSE,
                     n_chains = 3, burn = 2000, n_iter = 5000) {
  
  # jdat <- as.list(dat[, .(Y, A, B, AB, A_c, B_c, AB_c)])
  if (centered) {
    jdat <- as.list(dat[, .(Y, A_c, B_c, AB_c)])
    vars <- c("alpha", "gamma_a", "gamma_b", "gamma_ab")
  } else {
    jdat <- as.list(dat[, .(Y, A, B, AB)])
    vars <- c("alpha", "beta_a", "beta_b", "beta_ab")
  }
  jdat$N <- nrow(dat)
  
  mod <- jags.model(
    textConnection(model_string),
    data = jdat,
    n.chains = n_chains,
    quiet = TRUE
  )
  
  update(mod, burn, progress.bar = "none")
  
  samp <- coda.samples(
    mod,
    variable.names = vars,
    n.iter = n_iter,
    progress.bar = "none"
  )
  
  samp
}

samp_01 <- fit_jags(dd, model_01, centered = FALSE)
samp_c  <- fit_jags(dd, model_c, centered = TRUE)

diag_tbl <- function(samp, model_name) {
  post <- as_draws_df(samp)
  summ <- summarise_draws(post)
  out <- as.data.table(summ)
  out[, model := model_name]
  out[]
}

diag_01 <- diag_tbl(samp_01, "0/1-coded")
diag_c  <- diag_tbl(samp_c, "centered")

get_lor_summary <- function(samp, model_name) {
  dt <- as.data.table(as_draws_df(samp))
  
  if (model_name == "0/1-coded") {
    dt[, lOR_A := beta_a]
    dt[, lOR_B := beta_b]
    dt[, lOR_AB := beta_a + beta_b + beta_ab]
  } else {
    dt[, lOR_A := gamma_a - 0.5 * gamma_ab]
    dt[, lOR_B := gamma_b - 0.5 * gamma_ab]
    dt[, lOR_AB := gamma_a + gamma_b]
  }
  
  dt[, .(
    mean_A = mean(lOR_A),
    mean_B = mean(lOR_B),
    mean_AB = mean(lOR_AB),
    sd_A = sd(lOR_A),
    sd_B = sd(lOR_B),
    sd_AB = sd(lOR_AB)
  )]
}

lor_01 <- get_lor_summary(samp_01, "0/1-coded")
lor_c  <- get_lor_summary(samp_c,  "centered")

one_run <- function(
  n = 2000,
  truth = c(alpha = -0.8, beta_a = 0.5, beta_b = 0.9, beta_ab = -0.3),
  n_chains = 3,
  burn = 1000,
  n_iter = 3000
) {
  
  dd <- s_gen(
    n = n,
    alpha = truth["alpha"],
    beta_a = truth["beta_a"],
    beta_b = truth["beta_b"],
    beta_ab = truth["beta_ab"]
  )
  
  samp_01 <- fit_jags(
    dd, model_01, centered = FALSE, 
    n_chains = n_chains, burn = burn, n_iter = n_iter
  )
  
  samp_c <- fit_jags(
    dd, model_c, centered = TRUE, 
    n_chains = n_chains, burn = burn, n_iter = n_iter)
  
  get_metrics <- function(samp, model_name) {
    post <- as_draws_df(samp)
    summ <- as.data.table(summarise_draws(post))
    summ[, model := model_name]
    summ[]
  }
  
  out <- rbindlist(list(
    get_metrics(samp_01, "0/1-coded"),
    get_metrics(samp_c,  "centered")
  ))
  
  out[]
}

nsim <- 500

sim_res <- rbindlist(mclapply(seq_len(nsim), function(i) {
  out <- one_run()
  out[, sim := i]
  out[]
}, mc.cores = 5))

I’m currently on a quest to better know and understand treesitter-based tooling for R. To make it short, treesitter is a tool for parsing code, for instance recognizing what is a function, an argument, a logical in a string of code. With tools built upon treesitter you can search, reformat, lint and fix, etc. your code. Exciting stuff, running locally and deterministically on your machine.

Speaking of “etc.”, Etienne Bacher helpfully suggested I also look at treesitter-based tooling for other languages to see what’s still missing in our ecosystem. This is how I stumbled upon difftastic by Wilfred Hughes, “a structural diff tool that understands syntax”. This means that difftastic doesn’t only compare line or “words” but actual syntax by looking at lines around the lines that changed (by default, 3), Even better, it understands R out of the box¹.

Many thanks to Etienne Bacher not only for making me discover difftastic but also for useful feedback on this post!

Installing difftastic

To install difftastic I downloaded a binary file for my system from the releases of the GitHub repository, as documented in the manual.

difftastic on two files

You can run difftastic on two files, a bit like you would use the waldo R package on two objects.

Let’s compare:

a <- gsub("bad", "good", x)

to

a <- stringr::str_replace(x, "bad", "good")

respectedly saved in old.R and new.R. The CLI is called difft not difftastic. I use the “inline” display rather than the two columns default in order to save horizontal space.

difft old.R new.R --display inline

We’d get to this nice looking diff:

The parentheses and "bad" and "good" arguments are ignored.

We can also get the JSON version of this diff, which is an unstable feature which usage requires setting an environment variable:

export DFT_UNSTABLE=yes
difft old.R new.R --display json

This gets us

{"aligned_lines":[[0,0],[1,1]],"chunks":[[{"lhs":{"line_number":0,"changes":[{"start":5,"end":9,"content":"gsub","highlight":"normal"},{"start":23,"end":24,"content":",","highlight":"normal"},{"start":25,"end":26,"content":"x","highlight":"normal"}]},"rhs":{"line_number":0,"changes":[{"start":5,"end":12,"content":"stringr","highlight":"normal"},{"start":12,"end":14,"content":"::","highlight":"keyword"},{"start":14,"end":25,"content":"str_replace","highlight":"normal"},{"start":26,"end":27,"content":"x","highlight":"normal"},{"start":27,"end":28,"content":",","highlight":"normal"}]}}]],"language":"R","path":"content/post/2026-03-26-difftastic/new.R","status":"changed"}

Now, none of this isn’t very useful because I would never compare files in this way… I use version control!

difftastic with Git

We can set difftastic as the external diff tool for Git globally or for the current project.

For instance with the gert R package, to set it locally:

gert::git_config_set("diff.external", "difft")

If I want to use the inline display I’d set:

gert::git_config_set("diff.external", "difft --display inline")

Then git diff will by default use difftastic. Most interestingly for me, git show --ext-diff will use difftastic. I never use git diff directly but I do look at more or less recent commits a lot.

Say I am interested in the commit that removed roxygen2’s dependency on stringi, I’ll run:

git show 7a1dd39866699a2b0a034bb15244c07698a1e2e7 --ext-diff

and get:

This isn’t spectacular because this is a small diff, but I enjoy the highlighting of the parentheses of the removed nested call, and of the logical.

Cool features of difftastic

Building on two examples of the difftastic homepage…

Ignoring formatting changes

Since formatters can so helpfully apply your formatting preferences, reviewing formatting changes in a patch that’s about something else entirely is useless and annoying. Imagine having a function definition that fits on a single line, then adding one argument to it.

Going from

f <- function(myarg1 = foo, myarg2 = bar) {}

to

f <- function(
  myarg1 = foo,
  myarg2 = bar,
  myarg3 = baz
) {}

Because the definition is now longer than 80 characters, your formatter might switch the definition to be on multiple lines. But the actually interesting change is the addition of one argument.

Native Git diff² would show:

Git with difftastic would show:

The matching of delimiters is why I found the difftastic’s display of the roxygen2 commit more pleasing.

Matching delimiters in wrappers

The Git diff can look a bit ugly when you simply move code from one function to the other.

Say we go from

f <- function() {
  1 + 1
}

to

f <- function() {
  g()
}

g <- function() {
  1 + 1
}

Git diff would show:

Whereas Git with difftastic would show:

Will I use difftastic?

I really like the concept behind difftastic and the few Git commits I looked at with it rendered nicely. Now, what’s missing for me to use difftastic a lot is its integration with the tools where I actually use Git:

• Positron including the GitLens extension;
• GitHub Pull Request Files tab.

In any case, I’ll continue learning about tools based on treesitter, some of which like Air and Jarl I can already use directly from my IDE.

Dear rOpenSci friends, it’s time for our monthly news roundup! You can read this post on our blog. Now let’s dive into the activity at and around rOpenSci!

rOpenSci HQ

rOpenSci Dev Guide 1.0.0: Trilingual and Improved

rOpenSci Software Peer Review’s guidance is gathered in an online book that keeps improving! It is now available in English, Spanish and Portuguese. Read more in the release announcement

Champions Program Update

We are still going through the Champions selection process, and we’re excited to share that the new group of mentors has already been selected and is now actively reviewing Champions applications.

This cohort brings together a wonderful mix of returning Champions stepping into mentorship roles, mentors continuing their contributions, and new members joining the program. The 2026 mentors are Andrea Gómez Vargas, Pablo Paccioretti, Alber Hamersson Sánchez Ipia, Erick Isaac Navarro Delgado, Francisco Cardozo, Luis Verde Arregoitia, Monika Ávila Márquez, Guadalupe Pascal, Pao Corrales, and Elio Campitelli. Together, they represent a diverse and vibrant community across Colombia, Mexico, Argentina, Brazil, and Bolivia, with some currently based in Switzerland, Canada, the United States, and Australia. We’re very happy to see this growing, interconnected network supporting the next cohort of Champions.

R-Universe update

You can now download artifacts and log files from R-Universe without being logged in with a GitHub account, for example https://ropensci.r-universe.dev/opencv#checktable.

Software review and usage of AI tools

Authors submitting new software for peer review are now required to explain potential usage of generative AI tools in their package development. All submission templates now include a mandatory check-box:

- [ ] Generative AI tools were used to produce some of the material in this submission.
If so, please describe usage, and include links to any relevant aspects of your repository.

This is the start of our updates to accommodate generative AI tools in package development, as described in our recent blog post. The next phase will involve updates to our Dev Guide, explaining requirements and recommendations for authors, reviewers, and editors. All updates are intended to permit generative AI tools to be used in any useful way, while minimising the burden on those who volunteer their own time to keep our software peer review service running.

Software review bot updates

The ropensci-review-bot now delivers an initial report to all new software pre-submissions and submissions, identifying the five most similar packages from both all rOpenSci packages, and all CRAN packages. The matches are generated by our ropensci-review-tools/pkgmatch package (itself reviewed in this review issue). Matching is based on an “term frequency-inverse document frequency” algorithm, using inverse document frequencies from all rOpenSci and CRAN packages. Similar package reports can also be manually triggered (by editors only) with @ropensci-review-bot similar packages, like in this example for the pkgmatch package itself.

Coworking

Read all about coworking!

• Tuesday April 7th 2026, 9:00 Americas Pacific (16:00 UTC) “Getting to know the CSID Network” with Steffi LaZerte and cohosts Irene Ramos and Adamu Saleh Saidu.
  • Learn more about the CSID Network
  • Meet cohosts, Irene Ramos and Adamu Saleh Saidu, and learn more about the CSID Network and how you might get involved.
• Tuesday May 5th 2026, 9:00 Australia Western (01:00 UTC) “Code Review with rOpenSci” with Steffi LaZerte and cohost Liz Hare.
  • Explore resources for Code Review
  • Sign up to volunteer to do software peer-review at rOpenSci
  • Meet cohost, Liz Hare, and discuss resources for Code Review with rOpenSci.

And remember, you can always cowork independently on work related to R, work on packages that tend to be neglected, or work on what ever you need to get done!

Software

New packages

The following package recently became a part of our software suite:

• suwo, developed by Marcelo Araya-Salas together with Jorge Elizondo-Calvo and Alejandro Rico-Guevara: Streamline searching/downloading of nature media files (e.g. audios, photos) from online repositories. The package offers functions for obtaining media metadata from online repositories, downloading associated media files and updating data sets with new records. It has been reviewed by Eric R. Scott and Hugo Gruson.

Discover more packages, read more about Software Peer Review.

New versions

The following eleven packages have had an update since the last newsletter: cffr (v1.3.0), pkgmatch (v0.5.2), tarchetypes (0.14.1), rgbif (v3.8.5), saperlipopette (v0.1.1), gutenbergr (v0.5.0), trud (v0.2.1), naijR (v0.7.0), sasquatch (v0.1.3), lingtypology (v1.1.25), and rerddap (v1.2.3).

Post on dfms release: Releasing dfms 1.0: Fast and Feature-Rich Estimation of Dynamic Factor Models in R.

Software Peer Review

There are fifteen recently closed and active submissions and 5 submissions on hold. Issues are at different stages:

• One at ‘6/approved’:

  • suwo, Access Nature Media Repositories Through R. Submitted by Marcelo Araya-Salas.

• One at ‘5/awaiting-reviewer(s)-response’:

  • pkgmatch, Find R Packages Matching Either Descriptions or Other R Packages. Submitted by mark padgham.

• Two at ‘4/review(s)-in-awaiting-changes’:

  • logolink, An Interface for Running NetLogo Simulations. Submitted by Daniel Vartanian.

  • galamm, Generalized Additive Latent and Mixed Models. Submitted by Øystein Sørensen. (Stats).

• Six at ‘3/reviewer(s)-assigned’:

  • pvEBayes, Empirical Bayes Methods for Pharmacovigilance. Submitted by Yihao Tan. (Stats).

  • saperlipopette, Create Example Git Messes. Submitted by Maëlle Salmon.

  • ernest, A Toolkit for Nested Sampling. Submitted by Kyle Dewsnap. (Stats).

  • rcrisp, Automate the Delineation of Urban River Spaces. Submitted by Claudiu Forgaci. (Stats).

  • reviser, Tools for Studying Revision Properties in Real-Time Time Series Vintages. Submitted by Marc Burri.

  • priorsense, Prior Diagnostics and Sensitivity Analysis. Submitted by Noa Kallioinen. (Stats).

• Two at ‘2/seeking-reviewer(s)’:

  • nycOpenData, Convenient Access to NYC Open Data API Endpoints. Submitted by Christian Martinez.

  • RAMEN, RAMEN: Regional Association of Methylome variability with the Exposome and geNome. Submitted by Erick Navarro-Delgado.

• Three at ‘1/editor-checks’:

  • RAQSAPI, A Simple Interface to the US EPA Air Quality System Data Mart API. Submitted by mccroweyclinton-EPA.

  • fcmconfr, Fuzzy Cognitive Map Analysis in R. Submitted by benroston. (Stats).

  • coevolve, Fit Bayesian Generalized Dynamic Phylogenetic Models using Stan. Submitted by Scott Claessens. (Stats).

Find out more about Software Peer Review and how to get involved.

On the blog

Software Review

• Software Review in the Era of AI: What We Are Testing at rOpenSci by Mark Padgham, Noam Ross, Maëlle Salmon, Yanina Bellini Saibene, Mauro Lepore, Emily Riederer, Jouni Helske, and Francisco Rodriguez-Sanchez. rOpenSci is testing preliminary policies on the use of generative AI tools, with proposed updates to documentation and procedures for authors submitting software for review, for editors, and for reviewers.

• rOpenSci Dev Guide 1.0.0: Trilingual and Improved by Maëlle Salmon, Mark Padgham, and Noam Ross. Updates in version 1.0.0 of the online book ‘rOpenSci Packages: Development, Maintenance, and Peer Review’. Other languages: rOpenSci Dev Guide 1.0.0: Trilingüe y mejorada (es), Guia de desenvolvimento da rOpenSci 1.0.0: trilíngue e aprimorado (pt).

• Breaking Release of the patentsview R Package by Russ Allen and Chris Baker. Breaking Release of the patentsview R Package.

Calls for contributions

Calls for maintainers

If you’re interested in maintaining any of the R packages below, you might enjoy reading our blog post What Does It Mean to Maintain a Package?.

• NLMR, R package to simulate neutral landscape models. Issue for volunteering.

• landscapetools, R package for some of the less-glamorous tasks involved in landscape analysis. Issue for volunteering.

• hddtools, Tools to discover hydrological data, accessing catalogues and databases from various data providers. Issue for volunteering.

• qualtRics, download Qualtrics survey data. Issue for volunteering.

Calls for contributions

Refer to our help wanted page – before opening a PR, we recommend asking in the issue whether help is still needed.

Package development corner

Some useful tips for R package developers.

A new R core member!

The R Foundation announced that Heather Turner has joined the R Core Team!

How to browse the R mailing lists

The official mailing lists of the R project like R-package-devel are full of important and useful information. How to browse them, given that the default website is not easy to search? You can use the mail-archive website (thanks to Hugo Gruson for the reminder!) or a new project by James Balamuta: the R Mailing Lists Archive!

“Claude Code: Setting up ast-grep with R support”

Thanks to Mauro Lepore for sharing this blog post by Emil Hvitfeldt: “Claude Code: Setting up ast-grep with R support”. ast-grep is a tool for querying code by syntax rather than brittle regular expressions. The blog post describes how to add R support to this tool, and how to take advantage of it when using Claude.

On muffling messages from packages

A follow-up on our post “Please Shut Up! Verbosity Control in Packages”.

• With the {cli} R package you can change the default handler for messages. See the docs. It seems mostly used to muffle messages, e.g. in flir.
• Here’s how the usethis R package muffles gert message selectively.

Last words

Thanks for reading! If you want to get involved with rOpenSci, check out our Contributing Guide that can help direct you to the right place, whether you want to make code contributions, non-code contributions, or contribute in other ways like sharing use cases. You can also support our work through donations.

If you haven’t subscribed to our newsletter yet, you can do so via a form. Until it’s time for our next newsletter, you can keep in touch with us via our website and Mastodon account.

With the conflict in Iran causing worldwide disruption to energy markets, I have both a work and personal interest in energy supply in Pacific islands, which led me to this blog post. Here I look at just two aspects of energy: electricity generation, and household cooking. Nothing fancy here, just accessing some data and drawing a couple of plots.

Electricity generation

Here is the source of electricity for Pacific island countries, plus Australia and New Zealand, collated by Our World In Data from Energy Institute data that ultimately comes from government estimates:

There’s a pretty obvious story here: most of the Pacific is very dependent on “oil” (in the form of diesel) for generation of most of its electricity. There are some small steps towards renewables happening in recent years, but the vulnerability to a price or availability shock for diesel is pretty obvious.

Here’s the code for producing that, using the valuable owidapi R package to access the Our World in Data API.

#---------------Set up-----------------
library(owidapi)
library(tidyverse)
library(countrycode)
library(WDI)
library(jsonlite)
library(janitor)
library(httr2)


pic_codes <- 
  c(
    "ASM", "COK", "FSM", "FJI", "PYF", "GUM", "KIR", "MHL", "NRU", "NCL",
    "NIU", "MNP", "PLW", "PNG", "PCN", "WSM", "SLB", "TKL", "TON", "TUV",
    "VUT", "WLF", "AUS", "NZL"
  )
stopifnot(length(pic_codes) == 24)

# visual check we've got the right country codes for the Pacific:
countrycode::countrycode(pic_codes, origin = "iso3c", destination = "country.name.en")

#=======================electricity source===================

palette <- c(
  coal = "brown",
  gas = "magenta",
  oil = "red",
  kerosene = "red",
  electricity = "purple",
  solar = "yellow",
  wind = "steelblue",
  hydro = "darkblue",
  bioenergy = "lightgreen",
  charcoal = "grey",
  biomass = "darkgreen",
  'other renewables' = "darkgreen"
)

#-------------------electricity mix-----------------
elec_mix <- owid_get(
  chart_id = "share-elec-by-source",
  entities  = pic_codes
)

elec_data <- elec_mix |> 
   rename(country = entity_name) |> 
   select(-entity_id) |> 
   gather(variable, value, -country, -year) |>
   filter(value != 0) |> 
   filter(year > 2001) |> 
   mutate(variable = gsub("_share_of_electricity__pct", "", variable, fixed = TRUE),
          variable =gsub("_", " ", variable),
          variable =gsub(" excluding bioenergy", "", variable),
          variable = fct_drop(variable)) |> 
   mutate(variable = fct_relevel(variable, c("bioenergy", "hydro", "other renewables"), after = Inf)) |> 
   mutate(country = fct_relevel(country, c("Australia", "New Zealand"), after = Inf)) |> 
  group_by(country) |> 
  mutate(prop_pc= sum(value[variable %in% c("oil", "gas") & year == max(year)]) 
         / sum(value[year == max(year)])) |> 
  ungroup() |> 
  mutate(country = fct_reorder(country, prop_pc))

# Draw chart
elec_data |> 
  ggplot(aes(x = year, y = value, fill = variable)) +
 facet_wrap(~country, ncol = 5) +
  geom_col() +
  scale_fill_manual(values = palette) +
  scale_y_continuous(label = percent_format(scale = 1)) +
   labs(y = "Percentage of electricity",
        fill = "Source:",
        title = "Share of electricity by source",
        subtitle = "Countries shown in increasing order of vulnerability of electricity to a petrochemicals price or availability crisis.",
        x = "",
        caption = "Source: Ember (2026); Energy Institute - Statistical Review of World Energy (2025). Data processed by Our World In Data.") +
   theme(axis.text.x = element_text(angle = 45, hjust = 1))

Cooking fuel

OK, so electricity generation could be threatened by a lack of diesel. What about household cooking? This next chart draws on the definitive World Health Organization Household Energy Database, which models (based on what household survey data that is available) what households are using to cook:

Again, we see a lot of reliance on petrochemical products, particularly kerosene and liquid natural gas. The latter has been promoted as a relatively clean and healthy fuel to cook with compared to burning biomass (e.g. wood, coconuts, etc).

The larger Melanesian countries, with high rural populations, are those with the greatest use still of biomass for cooking. Most Pacific island countries do most of their cooking with oil or gas derived energy (remembering from the first chart, that ‘electricity’ often means diesel, ultimately).

Here’s the code to produce that chart. I used an LLM (I forget which) for the code to access the API itself, but I tested it and tweaked it to match my style, and the chart of course is all my own code.

#--------------------cooking-------------------
# The definitive source is the WHO  WHO Household Energy Database 
# which draws on various household surveys
# See https://www.who.int/data/gho/data/themes/air-pollution/cooking-fuel-and-technology-database-by-fuel-category

# next half dozen lines of code were supplied by Co-pilot and minimally
# tweaked by me for my style
indicator_code <- "PHE_HHAIR_PROP_POP_CATEGORY_FUELS"  # % by fuel type [3](https://millenniumindicators.un.org/wiki/spaces/SDGeHandbook/pages/35291272/Indicator+7.1.2)
url <- paste0("https://ghoapi.azureedge.net/api/", indicator_code)

resp <- request(url) |> 
  req_headers(`Accept` = "application/json")  |> 
  req_perform()

cooking_data <- fromJSON(resp_body_string(resp), flatten = TRUE)$value |>
  as_tibble() |> 
  clean_names()


pic_cooking_data <-cooking_data |> 
  filter(spatial_dim %in% pic_codes) |> 
  filter(dim1 == "RESIDENCEAREATYPE_TOTL") |> 
  mutate(fuel_type = tolower(gsub("HOUSEHOLDCOOKINGFUEL_FUEL_", "", dim2))) |> 
  mutate(year = as.numeric(time_dimension_value)) |> 
   select(value = numeric_value, 
          iso3_code = spatial_dim,
          value = numeric_value,
          year,
          fuel_type) |> 
  mutate(country = countrycode(iso3_code, origin = "iso3c", destination = "country.name.en"),
         country = gsub("Federated States", "Fed St", country)) |>
  group_by(country) |> 
  mutate(prop_gke = sum(value[fuel_type %in% c("gas", "kerosene", "electricity") & year == max(year)]) 
         / sum(value[year == max(year)])) |> 
  ungroup() |> 
  mutate(country = fct_reorder(country, prop_gke))

# Draw chart
pic_cooking_data |> 
  ggplot(aes(y = value, x = year, fill = fuel_type)) +
  facet_wrap(~country, ncol = 5) +
  # the numbers don't add up to 100 always, due to being modelled estimates
  #, not fully MECE, not counting dual fuels, etc. Good practice advice
  # is to not force them to add to 100%
  geom_area() +
  scale_fill_manual(values = palette) +
  scale_y_continuous(label = percent_format(scale = 1)) +
  labs(title = "Household primary fuel used for cooking",
       subtitle = "Estimates are modelled by WHO, and not adding up to 100% is a known limitation.
Countries shown in increasing order of vulnerability of cooking to a petrochemicals price or availability crisis.",
       x = "",
       fill = "Fuel type:",
       y ="Proportion of households",
       caption = "Source: WHO Household Energy Database")

That’s all, just a quick one today.

Dipping my toes into financial statements — income, balance sheet & cash flow. Still don’t fully get it, but slowly piecing together the story these numbers tell. Warren Buffett makes it look easy Baby steps!

Motivations

I’ve always wanted to learn financial statement, what it means, what it tells us, what Warren Buffett sees in them. Following the book Warren Buffett and the Interpretation of Financial Statements: The Search for the Company with a Durable Competitive Advantage and Datacamp: Analyzing Financial Statement in Python, I’ve made some notes for myself and also create the metrics functions, so that I can use and view them easily in the future. I’ll be honest, I still don’t fully understand it, but at least I can refer back to this as I look at these statements more frequently.

Disclaimer:

This is purely for educational purposes. This is not a financial advice, nor am I a financial advisor. This is a note to myself. If you find any mistakes or error, please let me know. Thanks! Also, there are a lot of information on each section, I won’t be covering all of them, just mostly the metrics from the book and also points I found interesting.

Objectives:

• The Skeleton of Financial Statements
• Let’s Take An Example
• Income Statement
  • Gross Profit Margin
  • Depreciation
  • Interest Payment to Operating Income
  • Income Before Tax
  • Income After Tax
  • Net Earnings
  • Per Share Earning
  • Operating Margin
  • Metrics
• Balance Sheet
  • Current Assets
  • Net Receivables To Gross Sale Ratio
  • The Current Ratio
  • Property, Plant, and Equipment
  • Short Term Debt
  • Long Term Debt
  • Retained Earnings
  • Treasury Stock
  • Return On Shareholders’ Equity
  • Metrics
• Cash Flow
  • Operating Income
  • Capital Expenditure
  • Stock Buyback
  • Metrics
• Combine All Metrics
• Let’s Look At Another Examples
• Oppotunities For Improvement
• Lessons Learnt

The Skeleton of Financial Statements

A financial statement is a formal record that shows a company’s financial activities and position, typically consisting of three core components: the income statement (which shows revenues earned and expenses incurred to calculate profit or loss over a period), the balance sheet (which presents what the company owns as assets, what it owes as liabilities, and the difference between them as equity at a specific point in time), and the cash flow statement (which tracks the actual movement of cash in and out of the business through operating activities, investing activities, and financing activities). The income statement reveals profitability, and the cash flow statement shows liquidity and how money actually moves through the business.

If we were to think of a kid’s lemonaid shop, the income statement would show how much money the shop made from selling lemon tea and how much it spent on ingredients, paying Johnny hourly to sell (salary) to calculate the profit. The balance sheet would list the shop’s assets (like cash in the register, inventory of lemons, sugar, and any equipment) and liabilities (like loans or unpaid bills – money your parent you borrowed from to buy all of the above) to show the net worth of the business at a given moment. The cash flow statement would track the actual cash coming in from customers and going out for expenses, giving insight into whether the shop has enough liquidity to cover its day-to-day operations.

It sounds simple, in the big picture, but these are just the basic skeleton of financial statements. There are many nuances and details that we need to understand to really grasp the story that these statements are telling us.Each section has its own items and some of these items are good at forming different metrics to tell the story of how the lemonaid business is doing. Below is just a snapshot of Apple’s financial statement.

Income Statement

Balance Sheet

Cash Flow

Let’s Take An Example

Let’s go to Alpha Vantage and create a free api key and then pull Apple’s 10 year financial statement and go through as an exercise.

https://www.sec.gov/Archives/edgar/data/320193/000032019325000079/aapl-20250927.htm#i719388195b384d85a4e238ad88eba90a_181

library(httr)
library(jsonlite)
library(tidyverse)

api_key <- "your_api_key_here" # Use .Renviron to be safer like below combine_all code

## Create a function to pull data
get_data <- function(fx,ticker) {
  raw <- GET(paste0(
    "https://www.alphavantage.co/query?function=",fx,
    "&symbol=",ticker,"&apikey=", api_key
  )) %>%
    content(as = "text", encoding = "UTF-8") %>%
    fromJSON()
  
  df <- raw$annualReports |> 
    as_tibble() |> 
    mutate(across(-c(fiscalDateEnding, reportedCurrency), as.numeric)) |>
    mutate(fiscalDateEnding = as.Date(fiscalDateEnding)) |>
    arrange(fiscalDateEnding)

  return(df)
}

## financial statement
income <- get_data("INCOME_STATEMENT","AAPL")
balance <- get_data("BALANCE_SHEET","AAPL")
cashflow <- get_data("CASH_FLOW","AAPL")

Let’s visualize the income statement

Income Statement

Alright, here we go! Income statement tells the story of how the goods and services are doing in the market, how much it costs to produce and sell them, and how much profit is left after all expenses are accounted for. It is a dynamic statement that shows the flow of money over a period of time, typically a quarter or a year. It is like a movie that tells the story of the company’s operations and profitability. Below are some ratios and heuristics of which companies have a durable competitive advantage, according to Warren Buffett.

Gross Profit Margin

\(Gross Profit Margin = Gross Profit / Total Revenue\)

In the book, we are looking for consistent, as a general rule, above 40% to consider a company to have a durable competitive advantage. Gross profit here is basically Total Revenue - Cost of Revenue, does not take into account of R&D, admin etc. Gross profit margin is a measure of how much profit a company makes from its core operations, before accounting for other expenses. A higher gross profit margin indicates that the company has a strong competitive position in the market and is able to generate more profit from its sales.

So on Apple, you can see that in the sales section, there is product and service. I assume product is the hardware and service is like the cloud storage etc. The cost of revenue section has the same sections, that depicts how much cost to make these products/services. Again, remember this is all just about the goods, does not take into account of the R&D, office that manages these, administrative etc, I think.

Depreciation

Apparently this is a non-cash expense that reflects the reduction in value of a company’s assets over time. It is an accounting method used to allocate the cost of tangible assets (like machinery, equipment, buildings) and intangible assets (like patents, copyrights) over their useful lives. Depreciation allows companies to spread out the expense of an asset over several years, rather than recognizing the entire cost in the year it was purchased.

A quick ratio in the book is

\(Depreciation / Gross Profit\)

which should be low ~5-7% which indicates that the company is not heavily reliant on physical assets that may lose value over time. Warren stated that EBITDA (Earnings Before Depreciation, Taxes, and Amortization) is something Wall Street loves, but in Warren’s eyes, depreciation is a real expanse. Why is it that he said EBITDA is something Wall Street love? Because it shows a higher profit by excluding non-cash expenses like depreciation, which can make a company look more profitable than it actually is… interesting.

Interest Payment to Operating Income

\(Interest Payment / Operating Income\)

should be less than 15%. This ratio indicates that the company is not overly burdened by debt and has a healthy balance between its operating income and interest expenses. A lower ratio suggests that the company is generating sufficient operating income to cover its interest payments, which is a positive sign of financial stability.

Operating Income = Total Revenue - Cost of Revenue - Operating Expenses. And Operating Expenses = R&D + Selling, General and Administrative Expenses. Wow, so many terminologies and I still don’t fully understand them, but at least I can refer back to this when I look at these statements.

Income Before Tax

This is a value that Warren uses when he calculates the return that he is getting when he a whole business.

Income After Tax

This is a value is the truth test of a business. What is reported with SEC should reflect the pre-tax income reported on the income statement. Apparently some companies like to report higher values than the the truth? was what the book said

Net Earnings

\(Net Earning / Total Revenue\)

The heuristic for this is we should look for more than 20% which indicates that the company is able to generate a significant amount of profit from its total revenue, hence a long term competitive advantage.

Net Earnings = Operating Income - Interest Expense - Taxes. This is the bottom line of the income statement, which represents the company’s total profit after all expenses have been deducted from total revenue. It is a key indicator of a company’s profitability and financial performance. A higher net earnings figure indicates that the company is generating more profit from its operations, which can be a sign of a strong competitive position in the market.Interest Expense - Taxes`.

Per Share Earning

This is something Warren wants to see consistent increment over time.

This is calculated by: \(EPS = Net Income/Outstanding Share\)

Does that mean if net income is negative, we’ll have negative share?

Operating Margin

\(Operating Margin = Operating Income / Total Revenue\)

This is a measure of a company’s profitability that indicates how much profit it generates from its operations relative to its total revenue. A higher operating margin suggests that the company is more efficient at converting revenue into profit, which can be a sign of a strong competitive position in the market. Warren looks for more than 10% which indicates that the company has a durable competitive advantage. ?????

Metrics

income_metric <- function(df) {
  df_i <- df |>
    mutate(grossProfitMargin = grossProfit / totalRevenue,
           depreciationToGrossProfit = depreciationAndAmortization / grossProfit,
           interestExpenseToOperatingIncome = interestExpense / operatingIncome,
           netEarningMargin = netIncome / totalRevenue,
           operatingMargin = operatingIncome / totalRevenue)
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "grossProfitMargin", 0.4, "red", 
    "depreciationToGrossProfit", 0.07, "blue",
    "interestExpenseToOperatingIncome", 0.15, "blue",
    "netEarningMargin", 0.2, "red",
    "operatingMargin", 0.1, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_i |>
    select(fiscalDateEnding, grossProfitMargin:operatingMargin) |>
    pivot_longer(cols = c(grossProfitMargin:operatingMargin), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw()
  
  return(plot)
}

income_metric(income) 

Net margin

What is EBTDA? It stands for Earnings Before Depreciation, Taxes, and Amortization. It is a measure of a company’s operating performance that excludes non-operating expenses such as depreciation, taxes, and amortization. It is often used as a proxy for cash flow from operations, as it focuses on the core profitability of the business before accounting for non-cash expenses and tax obligations.

Balance Sheet

\(Asset = Liability + Shareholder Equity\)

Reminded me of Full Metal Alchemist famous phrase “tōka kōkan”, equivalent exchange.

Current Assets

Current asset is any asset that can be reasonably expected to be converted into cash within one year. This includes cash and cash equivalents, accounts receivable, inventory, and other short-term assets. Current assets are important because they provide insight into a company’s liquidity and ability to meet its short-term obligations.

Net Receivables To Gross Sale Ratio

\(Net Receivables / Gross Sale\)

If a company consistently has a lower ratio (? 5% or less), it may indicate that the company is efficient in collecting payments from customers and has a lower risk of bad debts. A higher ratio may suggest that the company is having difficulty collecting payments, which could lead to cash flow issues and potential losses from uncollected receivables.

The Current Ratio

\(Current Ratio = Total Assets / Total Liabilities\)

The current ratio is a liquidity ratio that measures a company’s ability to pay off its short-term liabilities with its short-term assets. A current ratio of 1 or higher is generally considered good, indicating that the company has enough assets to cover its liabilities. A current ratio below 1 may indicate that the company may have difficulty meeting its short-term obligations. This makes sense.

Property, Plant, and Equipment

This is the value of a company’s physical assets, such as land, buildings, machinery, and equipment. It is important to consider the value of PPE when evaluating a company’s financial health and potential for growth. A company with a significant amount of PPE may have a competitive advantage in its industry, as it may be able to produce goods or services more efficiently than its competitors. However, it is also important to consider the age and condition of the PPE, as well as any potential liabilities associated with it. In the title of this chapter, it says For Warren Not Having Them Is A Good Thing

Short Term Debt

From Warren’s perspective, when it comes to investing in financial institutions, he’s always shied away from companies who are bigger borrowers in short-term than long-term.

\(Short Term Debt/ Long Term Debt\)

Not really sure what the heuristic threshold is, but let’s use less than 1 as a good indicator?

Long Term Debt

As a general rule, a company with durable competitve edge, will have no or little long term debt to maintain their business operation. Let’s ensure this is not an uptrend on visualization

Retained Earnings

Retained earnings is the portion of a company’s net income that is retained and not distributed as dividends to shareholders. It represents the accumulated profits that a company has reinvested in its business over time. Retained earnings can be used for various purposes, such as funding research and development, expanding operations, paying off debt, or acquiring other companies. It is an important metric for investors to consider when evaluating a company’s financial health and growth potential, as it indicates how much profit the company has generated and how it has been utilized to support its long-term success.

\(Growth Rate = (Ending Retained Earnings/Beginning Retained Earnings)^{1/years}-1\)

Treasury Stock

Treasury stock refers to shares that a company has repurchased from its shareholders. These shares are held in the company’s treasury and are not considered outstanding shares. Treasury stock can be used for various purposes, such as to increase shareholder value, to have shares available for employee compensation plans, or to prevent hostile takeovers. When a company repurchases its own shares, it reduces the number of outstanding shares in the market, which can increase the value of the remaining shares and potentially boost earnings per share (EPS). However, it is important for investors to consider the reasons behind a company’s decision to buy back its own stock and how it may impact the company’s financial health and long-term growth prospects.

Return On Shareholders’ Equity

\(Return On Shareholders' Equity = Net Income / Shareholders' Equity\)

According to the book, high ROSE means come play How about stop and smell the rose? The book did not tell us exactly what the threshold is, but the companies of choice has about ~30-35%

Metrics

balance_metric <- function(df, income) {
  df_b <- df |>
    mutate(currentRatio = totalCurrentAssets / totalCurrentLiabilities,
           netReceivablesToGrossSale = currentNetReceivables / income$grossProfit,
           shortToLongTermDebt = shortTermDebt / longTermDebt,
           growthRateRetainedEarnings = (retainedEarnings / lag(retainedEarnings))^(1/n()) - 1,
           returnOnShareholdersEquity = income$netIncome / totalShareholderEquity)
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "currentRatio", 1, "red", 
    "netReceivablesToGrossSale", 0.05, "blue",
    "shortToLongTermDebt", 1, "blue",
    "growthRateRetainedEarnings", 0.05, "red",
    "returnOnShareholdersEquity", 0.3, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_b |>
    select(fiscalDateEnding, currentRatio:returnOnShareholdersEquity) |>
    pivot_longer(cols = c(currentRatio:returnOnShareholdersEquity), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw()
  
  return(plot)
}

balance_metric(balance, income)

Cash Flow

Operating Income

Cash flow from operating income starts with net income and then add back depreciation and amortization. This is because depreciation and amortization are non-cash expenses that reduce net income but do not actually involve a cash outflow. By adding them back, we can get a better picture of the actual cash generated by the company’s operations.

Capital Expenditure

Also known as CapEx, refers to the funds that a company uses to acquire, upgrade, and maintain physical assets such as property, buildings, technology, equipment, or machinery. CapEx is an important metric for investors to consider when evaluating a company’s financial health and growth potential, as it indicates how much the company is investing in its long-term success.

If we were to look at Apple’s cash flow statement, CapEx is payments for acquisition of proptery, plant, and equipment.

\(Capital Expenditure/Net Earning\)

And the heuristic is ~50% or less.

Stock Buyback

Stock buyback, also known as share repurchase, refers to a company’s practice of buying back its own shares from the open market. This can be done for various reasons, such as to increase shareholder value, to have shares available for employee compensation plans, or to prevent hostile takeovers. When a company repurchases its own shares, it reduces the number of outstanding shares in the market, which can increase the value of the remaining shares and potentially boost earnings per share (EPS).

On the Apple cash flow statement, stock buyback is listed as repurchase of common stocks. Others might use issuance of (retirement) stocks.

Metrics

cashflow_metric <- function(df) {
  df_c <- df |>
    mutate(operatingCashFlow = netIncome + depreciationDepletionAndAmortization,
           capitalExpenditureToNetEarning = capitalExpenditures / netIncome,
           stockBuybackToNetEarning = abs(proceedsFromRepurchaseOfEquity / netIncome))
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "operatingCashFlow", 0, "red", 
    "capitalExpenditureToNetEarning", 0.5, "blue",
    "stockBuybackToNetEarning", 0.5, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_c |>
    select(fiscalDateEnding, operatingCashFlow, capitalExpenditureToNetEarning, stockBuybackToNetEarning) |>
    pivot_longer(cols = c(operatingCashFlow, capitalExpenditureToNetEarning, stockBuybackToNetEarning), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw()
  
  return(plot)
}

cashflow_metric(cashflow)

Use key valuation metrics (P/E, EV/EBITDA, P/B, P/S, etc) to determine how cheap or expensive a stock is.

Combine All Metrics

library(httr)
library(jsonlite)
library(tidyverse)
library(ggpubr)

api_key <- Sys.getenv("avkey")

## Create a function to pull data
get_data <- function(fx,ticker,share=F) {
  raw <- GET(paste0(
    "https://www.alphavantage.co/query?function=",fx,
    "&symbol=",ticker,"&apikey=", api_key
  )) %>%
    content(as = "text", encoding = "UTF-8") %>%
    fromJSON()
  
  if (share==T) {
    df <- raw$annualEarnings |>
      select(fiscalDateEnding, reportedEPS) |>
      mutate(fiscalDateEnding = ymd(fiscalDateEnding))
  } else {
  df <- raw$annualReports |> 
    as_tibble() |>
    mutate(fiscalDateEnding = as.Date(fiscalDateEnding)) |>   
    mutate(across(where(is.character), as.numeric)) |>        
    arrange(fiscalDateEnding)
  }

  return(df)
}

income_metric <- function(df) {
  df_i <- df |>
    mutate(grossProfitMargin = grossProfit / totalRevenue,
           depreciationToGrossProfit = depreciationAndAmortization / grossProfit,
           interestExpenseToOperatingIncome = interestExpense / operatingIncome,
           netEarningMargin = netIncome / totalRevenue,
           operatingMargin = operatingIncome / totalRevenue)
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "grossProfitMargin", 0.4, "red", 
    "depreciationToGrossProfit", 0.07, "blue",
    "interestExpenseToOperatingIncome", 0.15, "blue",
    "netEarningMargin", 0.2, "red",
    "operatingMargin", 0.1, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_i |>
    select(fiscalDateEnding, grossProfitMargin:operatingMargin) |>
    pivot_longer(cols = c(grossProfitMargin:operatingMargin), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw() +
    ggtitle("Income Statement Metrics")
  
  return(plot)
}

share_metric <- function(df) {
    plot <- df |>
      mutate(reportedEPS = as.numeric(reportedEPS)) |>
      ggplot(aes(x = fiscalDateEnding, y = reportedEPS)) +
      geom_line() +
      geom_smooth() +
      theme_bw() +
      ggtitle("Share Metrics")
    return(plot)
}

balance_metric <- function(df, income) {
  df_b <- df |>
    mutate(currentRatio = totalCurrentAssets / totalCurrentLiabilities,
           netReceivablesToGrossSale = currentNetReceivables / income$grossProfit,
           shortToLongTermDebt = shortTermDebt / longTermDebt,
           growthRateRetainedEarnings = (retainedEarnings / lag(retainedEarnings))^(1/n()) - 1,
           returnOnShareholdersEquity = income$netIncome / totalShareholderEquity)
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "currentRatio", 1, "red", 
    "netReceivablesToGrossSale", 0.3, "blue",
    "shortToLongTermDebt", 1, "blue",
    "growthRateRetainedEarnings", 0.05, "red",
    "returnOnShareholdersEquity", 0.3, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_b |>
    select(fiscalDateEnding, currentRatio:returnOnShareholdersEquity) |>
    pivot_longer(cols = c(currentRatio:returnOnShareholdersEquity), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |>
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw() +
    ggtitle("Balance Sheet Metrics")
  
  return(plot)
}

cashflow_metric <- function(df, income) {
  df_c <- df |>
    mutate(operatingCashFlow = netIncome + depreciationDepletionAndAmortization,
           capitalExpenditureToNetEarning = capitalExpenditures / netIncome,
           stockBuybackToNetEarning = case_when(
             proceedsFromRepurchaseOfEquity < 0 ~ -proceedsFromRepurchaseOfEquity / income$netIncome,
             income$netIncome < 0 & proceedsFromRepurchaseOfEquity > 0 ~ NA_real_,
             income$netIncome < 0 & proceedsFromRepurchaseOfEquity < 0 ~-proceedsFromRepurchaseOfEquity / income$netIncome))
  
  hline <- tribble(
    ~param, ~hline_value, ~color2,
    "operatingCashFlow", 0, "red", 
    "capitalExpenditureToNetEarning", 0.5, "blue",
    "stockBuybackToNetEarning", 0.5, "red"
  ) |>
    mutate(
      param = factor(param),
      ymin = ifelse(color2 == "red", -Inf, hline_value),
      ymax = ifelse(color2 == "red", hline_value, Inf)
    )
  
  columns <- hline$param
  
  plot <- df_c |>
    select(fiscalDateEnding, operatingCashFlow, capitalExpenditureToNetEarning, stockBuybackToNetEarning) |>
    pivot_longer(cols = c(operatingCashFlow, capitalExpenditureToNetEarning, stockBuybackToNetEarning), 
                 names_to = "param", values_to = "values") |>
    mutate(param = factor(param, levels = columns)) |> 
    ggplot(aes(x = fiscalDateEnding, y = values)) +
    geom_rect(
      data = hline,
      aes(ymin = ymin, ymax = ymax, fill = color2),
      xmin = -Inf, xmax = Inf,
      alpha = 0.15,
      inherit.aes = FALSE
    ) +
    geom_line() +
    facet_wrap(. ~ param, scale = "free_y") +
    scale_fill_identity() +
    theme_bw() +
    ggtitle("Cash Flow Metrics")
  
  return(plot)
}

show_all <- function(name) {
  income <- get_data("INCOME_STATEMENT",name)
  share <- get_data("EARNINGS",name,share=T)
  balance <- get_data("BALANCE_SHEET",name)
  cashflow <- get_data("CASH_FLOW",name)
  plot1 <- income_metric(income) 
  plot2 <- share_metric(share)
  plot3 <- balance_metric(balance, income)
  plot4 <- cashflow_metric(cashflow, income)
  stacked_plot <- ggarrange(plot1,plot3,nrow=2)
  squished_plot <- ggarrange(plot4,plot2,ncol=2, widths = c(2,1))
  combineplot <- ggarrange(plotlist = list(stacked_plot,squished_plot), nrow = 2, heights = c(2,1))
  return(combineplot)
}

Let’s Look At Another Example

As a heuristics, we’ve made our dataviz where red is like lava, you want to stay above it. Blue is like the sky, you want to stay below it. The sweet zone is in between. When there is a loess curve, we are trying to see if the share is consistently increasing.

Apple

show_all("AAPL")

Looking the income metrics, Apple seems to be doing pretty good. Gross profit margin, net earning margin, operating income margin are all above thresholds. Depreciation to gross profit ratio is appropriately below the threshold, which means there is no recent purchase of property, plants, machines, etc where these are the ones that will add to depreciation. There is also a low to no interest expense, which is good! No debt? Next we move on to balance sheet metrics, current ratio is not great, liquidity is not great, net receivable is ?OK, I guess it make sense, if your products are popular and provide some sort of finance option, you might have some receivables. Short term debt is much lower than long term debt, which is good. The ROSE is smelling pretty good there too! In terms of cash flow metrics, it’s looking really good, it has high cash flowing in, low CapEx, and consistently buying its own share from 2020 onwards.

Microsoft

show_all("MSFT")

As for Microsoft, now that we’re a bit more familiar with the sections, let’s just string our read instead of separating them. Good consistent profit for all profit, operating and net income. Interestingly, depreciation is high along with CapEx, did they buy property, plants or machine? Ah, they built data center? Interest expense is low, which is good. Good liquidity factor with current ratio above threshold, net receivables is OK same as Apple (a competitor), short/long debt is good as well. ROSE is close to the threshold (worth watching). Cash is flowing, not much stock buyback, but overall EPS is consistently increasing. With its investment in data centers for Azure, next few years we should be seeing a persistent depreciation to gross profit ratio, along with CapEx, right? I read somewhere where you can’t place all depreciation in a single year.

NVIDIA

show_all("NVDA")

Now let’s take a look at NVIDIA. Wow, gross profit, operating income, and net earning margin all really good, better than Apple and Microsoft! Well, make sense with all these data centers, they need GPUs from NVIDIA. Look at that depreciation and CapEx, barely any! That’s great, less to maintain and the current plants they have are adequate to supply the demand. Low interest expense too, not much debt interest. Liquidity is great with high current ratio! Net receivables is good too, these big companies are paying NVIDIA back! ROSE is also very fragrant! Cash is flowing through the roof. Not much stock buyback. Also exponential EPS! You know, Warren did say that be mindful of companies with R&D cost.

Intel

show_all("INTC")

Wow, first glance with the profits, we’re seeing lava red. We then see high rise in depreciation and CapEx. Did they buy more property, plants, or machines? Intel announces plans for an investment of more than $28 billion for two new chip factories in Licking County, Ohio. And maybe Germany and Arizona too? Interesting. Notice that the interest expense towards the end of 2025 was in 2 digit but negative? It’s because the the operating income was negative and a high interest expense. I wonder if I should make this absolute as opposed to keeping the negative. Anyway, this means that not only was Intel in the red but also paying quite a bit of interest from the debt, I assume. The current ratio looks pretty good, same goes with the net receivables. Short term debt is also much lower compared to long term. ROSE not so good. As for the cash flow, when we add depreciation and demortization back to net income, it brought the red back up. Notice how the EPS dropped significantly lately, and notice that the buyback of share is missing? We code it where if it’s a positive value of buyback, which means it’s issuing stock instead of repurchasing, hence NA. Wow, what do you think? Will they be able to turn this around? Do data centers use Intel chips or AMD’s threadripper?

AMD

show_all("AMD")

Interesting income statement metric results. Good gross profit margin, but the net earning is barely. Interesting depreciation trend, what happened in 2024, where CapEx didn’t budge? Great liquidity, pretty high net receivables. Wow, very high short-to-long term debt in ?2021-2022. ROSE is essentially 0. With cash flowing good towards 2025, and also high buybacks 2022 ish. Finally, uptrending EPS! This is an odd one. Pasting my observation onto Claude and wow, these findings are due to Xilinx acquisition. That makes sense! Didn’t build new plants, old property, plants, machines have already been amortized, hence low depreciation. The acquisition is funded by debt, hence high short-to-long term debt ratio. Very interesting, indeed! So, it could be true that data centers prefer AMD over Intel given the good profit?

Google

show_all("GOOG")

Very good profit of all 3 metrics. After 2020, depreciation downtrended which is good and so did CapEx. Not much debt. Downtrending current ratio but still good. Very good net receivables! ROSE is coming up. Great cash flow. There is buyback of equity. Rising EPS. All great signs for Google! Being an all-service (no hardware?) company, this is quite good and healthy!

3M

show_all("MMM")

Now on to 3M, not a tech company, but let’s see if our workflow of looking at financial statement will help us tell a story. Profit is not as great as other tech companies we looked at, though gross profit is above threshold, but the net earnings is in the lava. Notice both net earnings and operating margins dipped quite low in 2024? And then a pretty high interest expense to operating income ratio in 2025? I wonder if because of losing money, they borrowed long term, hence short-to-long term debt is not too high? Current ratio is acceptable, but why is net receivables so high for 3M? Clients were not able to pay 3M? ROSE is pretty good in 2024 and 2025 even when profits weren’t. But why? Cash flow with high variance last 3 years. There is stock buyback in 2025. And even though EPS downtrended past 2 years but still remained high. This is a very interesting one as well. 3M should be a company where there’s competitive advantage because they make lots of daily usables which doesn’t need a whole lot of R&D. But why the anomaly in 2024 requiring debt? What do you think?

Disney

show_all("DIS")

Alright, what about Disney? Profit doesn’t seem as good as I expected. Both gross profit and net earning margins were below threshold, except for operating income margin. Interesting, in 2021, there is increase in depreciation, CapEx, and interest expense with low short-to-long debt ratio. It almost looked like they borrowed some money to purchase something new? Did they rebuild a park or something? In terms of liquidity, it’s on the lava zone in 2025 and seemed like it downtrended for the past 3 years as well. Net receivables is OK? Though I would think it should be lower? ROSE is red. Cash flow looks good, this is interesting because whatever that caused the depreciation when added back now they have cash. Also interesting to note that from 2020 to 2024 they were issuing stock rather than buying back. Their EPS appear to be quite volatile for the pat 3 years.

Boeing

show_all("BA")

Last but not least, let’s look at Boeing. Three profit margin metrics in the lava zone. Spike of interest expense in 2021 along with depreciation. Something physical was bought with borrowed money that’s paid long term. Current ratio is good. Essentially no net receivables. Consistent floating ROSE on the red sea. Operating cash flow in the red for a few years then positive in 2025. CapEx spiked in 2025, ?what was spent in capital where it doesn’t really depreciate? Buyback stock is interesting, missing several past few years, essentially issuing stocks when we have NA data. EPS downtrend to the negatives.

Opportunities For Improvement

• We need to recheck growth rate again, not sure if the code is correct
• perhaps build this into an R package, and develop further
• Learn about valuation metrics such as P/E, EV/EBITDA, P/B, P/S
• Not really sure if short-to-long term debt ratio of 1 is a good threshold, seems too high.
• include actual numbers on geom_label with ggrepel and reduce font size

Lessons learnt

• Depreciation and CapEx seem to have strong correlation, which makes sense.
• learnt Alpha Vantage API, quite straight forward
• learnt to look at financial statement and its metrics

If you like this article:

• please feel free to send me a comment or visit my other blogs
• please feel free to follow me on BlueSky, twitter, GitHub or Mastodon
• if you would like collaborate please feel free to contact me

Following the post on exact Shapley values for time series explainability, this post illustrates an example of how to use sensitivity analysis to explain time-series forecasts, based on the ahead::dynrmf model and external regressors. What is sensitivity analysis in this context? It’s about evaluating the impact of changes in the external regressors on the time-series forecast.

The post uses the ahead::dynrmf_sensi function to compute the sensitivities, and the ahead::plot_dynrmf_sensitivity function to plot the results.

First, install the package:

devtools::install_github("Techtonique/ahead")

Then, run the following code:

# devtools::install_github("Techtonique/ahead")
# install.packages(c("fpp2", "e1071", "patchwork"))

library(ahead)
library(fpp2)
library(patchwork)
library(e1071)

#' # Example 1: US Consumption vs Income
sensitivity_results_auto <- ahead::dynrmf_sensi(
y = fpp2::uschange[, "Consumption"],
xreg = fpp2::uschange[, "Income"],
h = 10
)

plot1 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_auto, 
                           title = "Sensitivity of Consumption to Income (Ridge)",
                           y_label = "Effect (ΔConsumption / ΔIncome)")

#' # Example 1: US Consumption vs Income
sensitivity_results_auto_svm <- ahead::dynrmf_sensi(
  y = fpp2::uschange[, "Consumption"],
  xreg = fpp2::uschange[, "Income"],
  h = 10, 
  fit_func = e1071::svm # additional parameter passed to ahead::dynrmf
)

plot2 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_auto_svm, 
                                        title = "Sensitivity of Consumption to Income (SVM)",
                                        y_label = "Effect (ΔConsumption / ΔIncome)")

 
# Example 2: TV Advertising vs Insurance Quotes
sensitivity_results_tv <- ahead::dynrmf_sensi(
 y = fpp2::insurance[, "Quotes"],
   xreg = fpp2::insurance[, "TV.advert"],
   h = 8
 )

plot3 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_tv,
                           title = "Sensitivity of Insurance Quotes to TV Advertising (Ridge)",
                           y_label = "Effect (ΔQuotes / ΔTV.advert)")

sensitivity_results_tv_svm <- ahead::dynrmf_sensi(
  y = fpp2::insurance[, "Quotes"],
  xreg = fpp2::insurance[, "TV.advert"],
  h = 8, 
  fit_func = e1071::svm # additional parameter passed to ahead::dynrmf
)

plot4 <- ahead::plot_dynrmf_sensitivity(sensitivity_results_tv_svm,
                                        title = "Sensitivity of Insurance Quotes to TV Advertising (SVM)",
                                        y_label = "Effect (ΔQuotes / ΔTV.advert)")

(plot1+plot2)

(plot3+plot4)

LLM-enhanced momentum investing combines traditional momentum signals with real-time news interpretation by large language models (LLMs). The idea is straightforward: stocks with strong past returns are candidates for momentum portfolios, but their inclusion and weight are refined by LLM-generated sentiment scores derived from firm-specific news. This hybrid approach improves risk-adjusted returns (Sharpe, Sortino) and is particularly effective in concentrated, high-conviction portfolios.

Key Parameters

1. Lookback Window (k)

• Definition: The number of past days of news considered for sentiment analysis.
• Role: Determines how much recent information the LLM uses to judge momentum continuation.
• Example: If k = 5, the model analyzes the last 5 business days of headlines and summaries for each stock.

2. Forecast Horizon (l)

• Definition: The period over which momentum continuation is predicted.
• Role: Sets the “target” for the LLM’s forecast — how far into the future the model should judge momentum persistence.
• Example: If l = 5, the LLM predicts whether momentum will continue for the next 5 trading days.
• Connection to Rebalancing: The forecast horizon typically aligns with the rebalancing cycle. For weekly rebalancing, the horizon is 5 days; for monthly rebalancing, it’s ~21 days.

3. Portfolio Size (m)

• Definition: The number of stocks selected after LLM scoring.
• Role: Controls how concentrated or diversified the portfolio is.
• Example: From the top 20 YTD performers, you might select the top 10 after sentiment scoring.

4. Rebalancing Frequency (T)

• Definition: How often the portfolio is updated with new signals.
• Role: Sets the rhythm of portfolio refresh — weekly, monthly, or quarterly.
• Example: Weekly rebalancing means recalculating momentum and sentiment scores every 5 trading days.

Concept

The strategy begins with a classic momentum screen: select the top 20 S&P 500 companies by year-to-date (YTD) performance. Instead of stopping there, the approach integrates large language model (LLM) sentiment analysis of firm-specific news. By analyzing the last 5 business days of headlines and summaries, the LLM produces a score indicating whether momentum is likely to continue.

These scores are then used to re-weight the portfolio, tilting allocations toward companies with stronger news sentiment. Finally, the portfolio is narrowed to the top 10 conviction stocks.

Selected Parameters

• Lookback Window: 5 days of firm-specific news.
• Rebalancing Frequency: Weekly updates of the portfolio.
• Forecast Horizon: 5 trading days (aligned with the rebalancing cycle).

This setup ensures that the LLM is asked to judge whether momentum will persist until the next rebalance, making the signals both short-term and actionable.

What the Code Does Step by Step

1. Fetching the Data
   • The R script first pulls all S&P 500 tickers.
   • It calculates YTD returns for each stock.
   • The top 20 stocks by performance are selected as momentum candidates.
2. News Sentiment Analysis with LLM
   • For each of these 20 stocks, the code queries Bing News for recent headlines with Azure AI Services.
   • The last 5 business days of news are collected.
   • These headlines and summaries are sent to a Microsoft Foundry-hosted LLM.
   • The LLM outputs a score (0–1) indicating whether sentiment supports momentum continuation or signals reversal.
3. Portfolio Tilting
   • LLM scores are normalized to [-1, +1].
   • Baseline equal weights are tilted according to these scores.
   • The top 10 stocks by adjusted weight form the final portfolio

4. Visualization

• A styled table is created using the gt package.
• Adjusted weights are color-coded (red–green gradient).
• The final portfolio is saved as an image (top10.png).

Strategic Insight

This approach mirrors the methodology in the Swiss Finance Institute paper:

• Momentum ranking provides the baseline.
• LLM sentiment scoring refines stock selection and weighting.
• Portfolio tilting integrates qualitative news signals into quantitative allocation.

library(httr)
library(jsonlite)
library(tidyquant)
library(tidyverse)
library(lubridate)
library(gt)
library(gtExtras)
library(scales)
library(showtext)
library(webshot2)

# 1. Environment & Auth Setup
sysfonts::font_add_google("Roboto Slab", "roboto_slab")
showtext_auto()

# Azure & Bing Credentials
bing_key           <- "<your-bing-key>"
bing_endpoint      <- "<your-bing-endpoint>" 
azure_llm_key      <- "<your-llm-key>"
azure_llm_endpoint <- "<your-llm-endpoint>"

# R Part (first): S&P 500 Screening with Monthly Returns
sp500_tickers <- 
  tq_index("SP500") %>% 
  select(symbol, company)

# Calculate YTD change
momentum_df <- 
  sp500_tickers %>%
  tq_get(get = "stock.prices", from = floor_date(today(), "year")) %>%
  group_by(symbol) %>%
  arrange(date) %>% # Ensure chronological order for first/last functions
  summarize(
    total_return = (last(adjusted) / first(adjusted)) - 1, 
    .groups = "drop"
  ) %>%
  inner_join(sp500_tickers, by = "symbol") %>%
  slice_max(total_return, n = 20) %>%
  select(symbol, company)

# R Part (second): News Search and LLM Analysis
analyze_momentum_continuation <- function(ticker, company_name) {
  
  # Construct the Bing News Search query: ticker + ' stock'
  query_str <- paste0(ticker, " stock")
  news_url <- paste0(bing_endpoint, "v7.0/news/search")
  
  # Call Bing News Search API
  news_res <- GET(news_url, add_headers(`Ocp-Apim-Subscription-Key` = bing_key), 
                  query = list(q = query_str, count = 5, freshness = "Day"))
  
  # Short pause to throttle calls and avoid rate limit errors
  Sys.sleep(1)
  
  news_text <- ""
  if (status_code(news_res) == 200) {
    content <- fromJSON(content(news_res, "text", encoding = "UTF-8"))
    if (length(content$value) > 0) {
      news_text <- paste(content$value$name, content$value$description, collapse = " | ")
    }
  }
  
  # Construct the LLM payload for Azure AI Foundry
  prompt_payload <- list(
    messages = list(
      list(role = "system", content = "You are an LLM Enhanced Momentum Investing Agent."),
      list(role = "user", content = paste0(
        "Headlines + Summaries for ", company_name, " (", ticker, "): ", news_text,
        "\nPerform sentiment analysis based on the last 5 days of news (lookback=5, horizon=5). ",
        "Infer whether sentiment supports momentum continuation or signals reversal. ",
        "Return a JSON object with two fields: 'subsector' (string) and 'llm_score' (string: probability 0-1)."))
    ),
    temperature = 0.1
  )
  
  llm_res <- POST(url = azure_llm_endpoint, 
                  add_headers(`api-key` = azure_llm_key, `Content-Type` = "application/json"),
                  body = prompt_payload, encode = "json")
  
  if (status_code(llm_res) == 200) {
    llm_out <- fromJSON(content(llm_res, "text", encoding = "UTF-8"))
    # Parse JSON response without using regex
    llm_json_data <- fromJSON(llm_out$choices$message$content)
    return(as.data.frame(llm_json_data))
  } else {
    return(data.frame(subsector = "N/A", llm_score = "0.5"))
  }
}

# Execute Analysis: Merge top 20 tickers with LLM scores
news_scores_df <- momentum_df %>%
  mutate(analysis = map2(symbol, company, analyze_momentum_continuation)) %>%
  unnest(analysis) %>%
  mutate(llm_score = as.numeric(llm_score))

# R Part (final): Portfolio Tilting and Visualization
# Normalize scores to [-1, +1] and tilt weights
tilted_portfolio <- news_scores_df %>%
  mutate(
    norm_score = rescale(llm_score, to = c(-1, 1), from = c(0, 1)),
    base_weight = 1 / n(),
    adj_weight = base_weight * (1 + norm_score)
  ) %>%
  mutate(adj_weight = adj_weight / sum(adj_weight)) %>%
  slice_max(adj_weight, n = 10)

# Create gt visualization using original column names
final_table <- tilted_portfolio %>%
  select(company, subsector, adj_weight) %>%
  gt() %>%
  tab_header(title = "Top 10 Tilted S&P 500 Momentum Stocks") %>%
  # Use cols_label for human-readable labels without renaming underlying columns
  cols_label(
    company = "Company", 
    subsector = "Subsector", 
    adj_weight = "Adjusted Weight "
  ) %>%
  # Apply color intensity with scales::col_numeric
  data_color(
    columns = adj_weight, 
    colors = col_numeric(palette = c("red", "green"), domain = NULL)
  ) %>%
  fmt_percent(
    columns = contains("adj_weight"), 
    decimals = 2,
    locale = "en" 
  ) %>% 
  cols_align(align = "center") %>%
  opt_table_font(font = google_font("Roboto Slab"))

# Save the visualization as top10.png using webshot
gtsave(final_table, "top10.png")

Final Observation

Looking at the resulting portfolio, one striking feature is the dominance of energy and petroleum companies. Firms such as ConocoPhillips, EOG Resources, ExxonMobil, Occidental Petroleum, Marathon Petroleum, Valero, Phillips 66, Chevron, and Baker Hughes all appear prominently.

This heavy tilt toward energy is not random—it reflects how geopolitical tensions (Iran and U.S.–Israel war dynamics) have amplified the importance of oil and gas in global markets. News sentiment around these companies has been strongly supportive of continued momentum, pushing them into the top allocation slots.

Over the course of my career as a Data Scientist, I’ve worked on projects ranging from simple code reviews, to large application builds. For the most part, I have used R to do this.

If you’re getting into coding or data science, one question you’re probably asking yourself is “Which language should I learn?”

This blog aims to show you why R might be a good decision.

R was built for data (not just programming)

Unlike general purpose languages (such as Python), R was designed specifically for statistics and data analysis.

That means:

• Built in statistical tools
• Powerful visualisation capabilities
• Research level methods available immediately

With packages like the tidyverse, you can clean, analyse, and visualise data with surprisingly little code.

High demand in analytics, research, and healthcare

R is especially popular in many sectors such as:

• Healthcare & biostats
• Academic research
• Government departments
• Finance & risk modeling
• Pharmaceutical companies

Here are some examples of R in production use:

• The {bbplot} R package. Yes, the BBC use R to create graphics for their website!
• Health and wellbeing profiling app for the NHS
• During the Covid-19 pandemic, the financial times had a stats tracker in which the graphs were built with R.

Knowing some R will give you a competitive edge if you’re looking at working within these sectors.

Open source with the backing of Posit

R is open source. This means that:

• It’s free, and always will be!
• Anyone can view the source code the makes up R, there are.
• Each R package (a folder containing code) has to live on GitHub.com, for everyone to see.
• It has a large community of contributors. There are great forums to get help such as Stack Overflow, Posit Community and the R weekly newsletter and tonnes more.
• There are thousands more available functionalities compared to paid softwares such as SPSS, SAS or Excel.

Posit, who maintain the free to use RStudio and Positron IDEs (integrated development environment), have many full time staff working solely on maintaining and creating new functionality within R. This means we get:

• Defined accountability
• Predictable release cycles
• Bugs can be solved quicker

Incredible data visualisation possibilities

Being able to communicate your findings with stakeholders is very important in data science, and one of R’s biggest strengths is visualisation and reporting.

With the {ggplot2} package, you can create publication ready charts with very little code. The R Graph Gallery has some amazing examples of what is possible with {ggplot2}.

With the {quarto} and {shiny} packages, you are able to build reproducible reports and interactive dashboards. All this without needing to know any HTML, CSS or JavaScript.

Beginner friendly learning curve

This is very much my own opinion. Compared to other languages, I think R is fairly intuitive and feels rewarding much earlier on in the journey. It also has (in my opinion), the most beginner friendly programme to code in, called RStudio.

Most people attend only two days worth of training with Jumping Rivers, and say they feel ready to start tackling their own data problems.

So… is R worth learning in 2026?

I think so. If you want pure software engineering or large-scale production systems, you may need Python. But for becoming a strong data thinker, and giving you an edge in your analysis, R is one of the best starting points.

For updates and revisions to this article, see the original post

x

Two hundred and fifty-five of the new packages submitted to CRAN in February were still there in mid-March. Here are my Top 40 picks in seventeen categories: Artificial Intelligence, Biology, Buddhism, Climate Science, Computational Methods, Data, Ecology, Epidemiology, Genomics, Machine Learning, Medical Applications, Physics, Statistics, Surveys, Time Series, Utilities, and Visualization.

In computational linguistics, word meanings are shaped by their contexts. As the British linguist John Rupert Firth put it in 1957, ‘You shall know a word by the company it keeps’ (see Brunila & LaViolette, 2022, for a re-examination of the intellectual history). It sounds almost like life advice, but Firth meant something technical: words that habitually appear alongside each other tend to share semantic territory. The adjective ‘good’, for instance, is far more likely to appear near ‘kind’, ‘genuine’, ‘fair’ and ‘quality’ than near ‘broken’ or ‘fraud’ – and a model that tracks those neighbours can learn what ‘good’ means without ever being told. The principle extends to polysemy: ‘bank’ means something entirely different in the company of ‘river’ and ‘fishing rod’ than in the company of ‘overdraft’ and ‘mortgage’. Context is everything.

This deceptively simple insight is the bedrock on which generative AI was built. The earliest computational implementations of Firth’s principle – distributional semantic models such as Latent Semantic Analysis (LSA; Landauer & Dumais, 1997) and the Hyperspace Analogue to Language (Lund & Burgess, 1996) – were modest by today’s standards: a matrix of word co-occurrence counts, a few hundred latent dimensions and a vocabulary of perhaps tens of thousands of words. Yet even these pocket-sized models captured real-world structure with startling fidelity. Louwerse and Zwaan (2009) showed that the frequency with which city names co-occur in English text predicts their actual geographical distances: cities close together on a map tend to be mentioned together more often, and an LSA model trained on text alone can reconstruct approximate maps of the United States without ever seeing one. Louwerse (2011) extended this further, showing that text statistics encode not just geography but sensory properties, emotional associations and conceptual relationships across a wide range of domains. Indeed, distributional language statistics may track some sensorimotor properties of concepts (Bernabeu, 2022; Louwerse & Connell, 2011; cf. Xu et al., 2025), especially after fine-tuning on human sensorimotor ratings (Wu et al., 2026). In short, language does not merely label the world – it encodes its structure, and even a simple co-occurrence model can read that encoding back.

We can see this for ourselves. The R code included below (click ‘Expand’ to view it) applies LSA – one of the simplest distributional models – to three text collections, projects the resulting word vectors into two dimensions via PCA (principal component analysis) and plots them. In brief, LSA builds a term-document matrix (a large table recording how often each word appears in each document), weights it with TF-IDF (term frequency–inverse document frequency, which highlights words distinctive to particular documents rather than ubiquitous everywhere) and then compresses it via truncated SVD (singular value decomposition, a form of dimensionality reduction). Each corpus is split into two groups: the most distinctive words per group (selected by the difference in mean TF-IDF weight between groups) are plotted in the group’s colour, while the most frequent shared words appear in purple. Words that co-occur in similar contexts cluster together; words from different domains drift apart.

PCA works by finding new axes – principal components – that capture the maximum variance in the data. Each word receives a loading on each component: a number ranging from −1 to +1 that indicates how strongly that word contributes to that axis of variation (a gentle introduction to PCA in R is available in an earlier post on this blog). High absolute loadings on a component mean that the word is a strong marker of the distinction that component captures.

How are the thematic groups decided? The code computes the mean TF-IDF weight of every word in each group of documents and then takes the difference. Words whose weight is much higher in group A than in group B are classified as distinctive to A, and vice versa. The top 15 words at each extreme become the coloured labels in the plot, while the most frequent words that do not belong to either extreme are labelled ‘Shared’. The grouping is therefore entirely data-driven: no human decides which words are ‘finance’ or ‘energy’ – the corpus statistics do. Above each plot, a table shows the mean loading of each thematic group on the first two principal components, with the highest positive loading per group highlighted in bold. A high absolute loading tells us that a given group of words is strongly aligned with that component – in other words, that the component captures precisely the distinction between those groups. When one group loads heavily on PC1 while another does not, the first principal component is essentially the axis that separates them.

pkgs <- c("LSAfun", "tm", "ggplot2", "plotly")
invisible(lapply(pkgs, function(p)
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p)))
library(LSAfun)
library(tm)
library(ggplot2)
library(plotly)

# --- Reusable helper: LSA + PCA plot ------------------------------------
# Builds a TF-IDF term-document matrix, computes a truncated SVD,
# selects the most distinctive and most shared words and projects them
# to 2D via PCA *on the selected words only* for maximum spread.

lsa_pipeline <- function(doc_list, labels, grp_a, grp_b,
                         lab_a, lab_b, colour_a, colour_b,
                         top_n = 15, n_shared = 10,
                         k = 20, min_docs = 4) {
  corp <- VCorpus(VectorSource(doc_list))
  corp <- tm_map(corp, content_transformer(tolower))
  corp <- tm_map(corp, removePunctuation)
  corp <- tm_map(corp, removeNumbers)
  corp <- tm_map(corp, removeWords, stopwords("en"))
  corp <- tm_map(corp, stripWhitespace)
  tdm  <- as.matrix(TermDocumentMatrix(corp,
             control = list(weighting = weightTfIdf,
                            bounds = list(global = c(min_docs, Inf)))))
  k_use <- min(as.integer(k), nrow(tdm) - 1L, ncol(tdm) - 1L)
  sv    <- svd(tdm, nu = k_use, nv = k_use)
  wlsa  <- sv$u %*% diag(sv$d[1:k_use])
  rownames(wlsa) <- rownames(tdm)
  idx_a  <- which(labels == grp_a)
  idx_b  <- which(labels == grp_b)
  mean_a <- rowMeans(tdm[, idx_a, drop = FALSE])
  mean_b <- rowMeans(tdm[, idx_b, drop = FALSE])
  total  <- mean_a + mean_b
  spec   <- mean_a - mean_b           # positive = distinctive to A
  top_a  <- names(sort(spec, decreasing = TRUE))[1:top_n]
  top_b  <- names(sort(spec, decreasing = FALSE))[1:top_n]
  shared_pool <- setdiff(names(sort(total, decreasing = TRUE)),
                         c(top_a, top_b))
  shared <- head(shared_pool, n_shared)
  hl   <- unique(c(top_a, top_b, shared))
  hl   <- hl[hl %in% rownames(wlsa)]
  # PCA on the selected words only, for better spatial spread
  wlsa_hl <- wlsa[hl, , drop = FALSE]
  pca  <- prcomp(wlsa_hl, scale. = FALSE)
  cd   <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2],
                     word = rownames(wlsa_hl))
  cd$topic <- ifelse(cd$word %in% top_a & !cd$word %in% top_b, lab_a,
              ifelse(cd$word %in% top_b & !cd$word %in% top_a, lab_b,
                     "Shared"))
  p <- ggplot(cd, aes(PC1, PC2, colour = topic,
                      text = paste0(word, " (", topic, ")"))) +
    geom_point(size = 0, alpha = 0) +
    scale_colour_manual(values = setNames(c(colour_a, colour_b, "#7B2D8E"),
                                          c(lab_a, lab_b, "Shared")),
                        guide = guide_legend(override.aes = list(size = 3, alpha = 1))) +
    labs(x = "Principal Component 1", y = "Principal Component 2",
         colour = NULL) +
    theme_minimal(base_size = 12) +
    theme(legend.position = "bottom",
          legend.margin   = margin(t = -5),
          axis.title.x    = element_text(margin = margin(t = 12)),
          axis.title.y    = element_text(margin = margin(r = 12)),
          plot.margin     = margin(0, 0, 0, 0))

  # Map each word to its group colour for label text
  col_map <- setNames(c(colour_a, colour_b, "#7B2D8E"),
                      c(lab_a, lab_b, "Shared"))
  cd$label_col <- col_map[cd$topic]

  # Trim spatial outliers so the dense cluster is readable.
  # Words beyond the IQR fence are dropped from the plot (not from LSA).
  q1  <- quantile(cd$PC1, 0.25); q3 <- quantile(cd$PC1, 0.75)
  iqr <- q3 - q1; fence <- 2.5
  keep <- cd$PC1 >= (q1 - fence * iqr) & cd$PC1 <= (q3 + fence * iqr)
  q1y <- quantile(cd$PC2, 0.25); q3y <- quantile(cd$PC2, 0.75)
  iqry <- q3y - q1y
  keep <- keep & cd$PC2 >= (q1y - fence * iqry) & cd$PC2 <= (q3y + fence * iqry)
  cd <- cd[keep, , drop = FALSE]

  pp <- ggplotly(p, tooltip = "text")
  # Hide all ggplot traces from plot AND legend
  for (k in seq_along(pp$x$data)) {
    pp$x$data[[k]]$marker$size    <- 0.1
    pp$x$data[[k]]$marker$opacity <- 0
    pp$x$data[[k]]$showlegend <- FALSE
  }
  # Constrain axes to the data range (with a small pad)
  pad_x <- diff(range(cd$PC1)) * 0.06
  pad_y <- diff(range(cd$PC2)) * 0.06
  # Add text traces per group (toggleable via legend)
  legend_groups <- c(lab_a, lab_b, "Shared")
  legend_cols   <- c(colour_a, colour_b, "#7B2D8E")
  offscreen_x <- max(cd$PC1) + pad_x * 50
  offscreen_y <- max(cd$PC2) + pad_y * 50
  for (i in seq_along(legend_groups)) {
    grp <- legend_groups[i]
    grp_data <- cd[cd$topic == grp, , drop = FALSE]
    if (nrow(grp_data) == 0) next
    # Text trace at actual positions (no legend entry)
    pp <- pp %>% add_trace(
      x = grp_data$PC1, y = grp_data$PC2,
      type = "scatter", mode = "text",
      text = grp_data$word,
      textfont = list(size = 11, color = legend_cols[i]),
      name = grp, legendgroup = grp, showlegend = FALSE,
      hoverinfo = "text",
      hovertext = paste0(grp_data$word, " (", grp, ")"),
      inherit = FALSE
    )
    # Legend-only marker trace (off-screen, linked via legendgroup)
    pp <- pp %>% add_trace(
      x = offscreen_x, y = offscreen_y, type = "scatter", mode = "markers",
      marker = list(size = 12, color = legend_cols[i], opacity = 1,
                    symbol = "circle"),
      name = grp, legendgroup = grp, showlegend = TRUE,
      hoverinfo = "skip", inherit = FALSE
    )
  }
  pp <- pp %>% layout(
    legend = list(orientation = "h", x = 1, xanchor = "right",
                  y = -0.12, tracegroupgap = 4, itemwidth = 30,
                  itemsizing = "constant",
                  font = list(size = 12),
                  bordercolor = "#CCCCCC", borderwidth = 1,
                  bgcolor = "#FAFAFA",
                  xpad = 4, ypad = 10),
    xaxis = list(title = list(text = "Principal Component 1",
                              standoff = 8),
                 range = c(min(cd$PC1) - pad_x, max(cd$PC1) + pad_x)),
    yaxis = list(title = list(text = "Principal Component 2",
                              standoff = 8),
                 range = c(min(cd$PC2) - pad_y, max(cd$PC2) + pad_y)),
    margin = list(b = 60)
  )
  list(plot = pp, lsa = wlsa, tdm = tdm, pca = pca, words = cd)
}

# --- 1. Reuters newswire ------------------------------------------------
data(acq)
data(crude)

docs   <- c(lapply(acq, content), lapply(crude, content))
labels <- c(rep("acq", length(acq)), rep("crude", length(crude)))

res1 <- lsa_pipeline(docs, labels,
  grp_a = "acq", grp_b = "crude",
  lab_a = "Finance", lab_b = "Energy",
  colour_a = "#D55E00", colour_b = "#0072B2",
  min_docs = 4)

# Cosine similarities in the 20-dimensional LSA space
pairs <- list(
  c("oil", "barrel"), c("shares", "acquisition"),
  c("price", "barrel"), c("price", "shares"),
  c("shares", "oil"), c("acquisition", "barrel"))
pairs <- Filter(function(p) all(p %in% rownames(res1$lsa)), pairs)
sims  <- sapply(pairs, function(p)
  round(Cosine(p[1], p[2], tvectors = res1$lsa), 3))
names(sims) <- sapply(pairs, paste, collapse = " ~ ")
sims
#>         oil ~ barrel shares ~ acquisition       price ~ barrel 
#>                0.675                0.236                0.938 
#>       price ~ shares         shares ~ oil acquisition ~ barrel 
#>                0.129               -0.014               -0.043

if (!requireNamespace("sotu", quietly = TRUE)) install.packages("sotu")

sotu_texts  <- sotu::sotu_text
sotu_years  <- sotu::sotu_meta$year
sotu_labels <- ifelse(sotu_years < 1945, "Pre-1945", "Post-1945")

res2 <- lsa_pipeline(as.list(sotu_texts), sotu_labels,
  grp_a = "Pre-1945", grp_b = "Post-1945",
  lab_a = "Pre-1945", lab_b = "Post-1945",
  colour_a = "#E69F00", colour_b = "#009E73",
  min_docs = 5)

if (!requireNamespace("text2vec", quietly = TRUE)) install.packages("text2vec")

data("movie_review", package = "text2vec")
mv_labels <- ifelse(movie_review$sentiment == 1, "Positive", "Negative")

res3 <- lsa_pipeline(as.list(movie_review$review), mv_labels,
  grp_a = "Positive", grp_b = "Negative",
  lab_a = "Positive", lab_b = "Negative",
  colour_a = "#009E73", colour_b = "#D55E00",
  min_docs = 50)

# Core data wrangling and visualization
library(tidyverse)

# Bioinformatics packages
library(DESeq2)
library(pheatmap)
library(limma)
library(edgeR)

# Statistical learning and modeling
library(caret)
library(glmnet)
library(randomForest)

# Time-aware analysis
library(forecast)
library(tsibble)
library(fable)

# Annotation and interpretation
library(clusterProfiler)
library(org.Hs.eg.db)

# Helpful utilities
library(broom)
library(ggrepel)
library(pROC)

set.seed(123)

theme_set(
  theme_minimal(base_size = 13) +
    theme(
      plot.title = element_text(face = "bold"),
      axis.title = element_text(face = "bold"),
      panel.grid.minor = element_blank()
    )
)

# Read count matrix and sample metadata
counts <- read.csv("gene_counts.csv", row.names = 1, check.names = FALSE)
metadata <- read.csv("sample_metadata.csv", row.names = 1)

# Ensure samples align
metadata <- metadata[colnames(counts), , drop = FALSE]

# Inspect dimensions
dim(counts)
dim(metadata)

# Preview data
head(counts[, 1:6])
head(metadata)

# Basic integrity checks
stopifnot(all(colnames(counts) == rownames(metadata)))
sum(is.na(counts))
sum(is.na(metadata))

# Explore metadata structure
str(metadata)
table(metadata$condition)
table(metadata$batch)
table(metadata$timepoint)

# Total reads per sample
library_sizes <- colSums(counts)
sort(library_sizes)

# Visualize library size distribution
library_df <- tibble(
  sample = names(library_sizes),
  total_counts = library_sizes
)

ggplot(library_df, aes(x = reorder(sample, total_counts), y = total_counts)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Library Size per Sample",
    x = "Sample",
    y = "Total Counts"
  )

# Filter low-count genes
keep_genes <- rowSums(counts >= 10) >= 3
counts_filtered <- counts[keep_genes, ]

dim(counts_filtered)

# Optional: identify highly variable genes after transformation
log_counts <- log2(counts_filtered + 1)
gene_variance <- apply(log_counts, 1, var)

hv_genes <- names(sort(gene_variance, decreasing = TRUE))[1:500]
length(hv_genes)
head(hv_genes)

# Build DESeq2 object
dds <- DESeqDataSetFromMatrix(
  countData = counts_filtered,
  colData = metadata,
  design = ~ batch + condition
)

# Run differential expression pipeline
dds <- DESeq(dds)

# Extract normalized counts
norm_counts <- counts(dds, normalized = TRUE)

# Results for treatment vs control
res <- results(dds, contrast = c("condition", "treated", "control"))

# Order by adjusted p-value
res_ordered <- res[order(res$padj), ]
res_df <- as.data.frame(res_ordered) %>%
  rownames_to_column("gene")

head(res_df)

# Summary of results
summary(res)

# Significant genes
sig_res <- res_df %>%
  filter(!is.na(padj), padj < 0.05, abs(log2FoldChange) > 1)

nrow(sig_res)
head(sig_res)

# Variance stabilizing transformation
vsd <- vst(dds, blind = FALSE)
vsd_mat <- assay(vsd)

# Principal component analysis
pca_data <- plotPCA(vsd, intgroup = c("condition", "batch"), returnData = TRUE)
percent_var <- round(100 * attr(pca_data, "percentVar"))

ggplot(pca_data, aes(PC1, PC2, color = condition, shape = batch, label = name)) +
  geom_point(size = 4) +
  geom_text_repel(size = 3.5, max.overlaps = 20) +
  labs(
    title = "PCA of Variance-Stabilized Expression Data",
    x = paste0("PC1: ", percent_var[1], "% variance"),
    y = paste0("PC2: ", percent_var[2], "% variance")
  )

# Sample-to-sample distance heatmap
sample_dists <- dist(t(vsd_mat))
sample_dist_matrix <- as.matrix(sample_dists)

pheatmap(
  sample_dist_matrix,
  clustering_distance_rows = sample_dists,
  clustering_distance_cols = sample_dists,
  main = "Sample Distance Heatmap"
)

# Volcano plot
volcano_df <- res_df %>%
  mutate(
    significance = case_when(
      !is.na(padj) & padj < 0.05 & log2FoldChange > 1  ~ "Upregulated",
      !is.na(padj) & padj < 0.05 & log2FoldChange < -1 ~ "Downregulated",
      TRUE ~ "Not significant"
    ),
    neg_log10_padj = -log10(padj)
  )

ggplot(volcano_df, aes(x = log2FoldChange, y = neg_log10_padj, color = significance)) +
  geom_point(alpha = 0.75) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  labs(
    title = "Volcano Plot of Differential Expression",
    x = "Log2 Fold Change",
    y = "-Log10 Adjusted P-value"
  )

# Heatmap of top significant genes
top_genes <- sig_res %>%
  slice_min(order_by = padj, n = 30) %>%
  pull(gene)

heatmap_mat <- vsd_mat[top_genes, ]
heatmap_mat_scaled <- t(scale(t(heatmap_mat)))

annotation_col <- metadata %>%
  select(condition, batch)

pheatmap(
  heatmap_mat_scaled,
  annotation_col = annotation_col,
  show_rownames = TRUE,
  show_colnames = TRUE,
  clustering_method = "complete",
  main = "Top Differentially Expressed Genes"
)

# Example: prepare a classification dataset using selected genes
selected_genes <- sig_res %>%
  slice_min(order_by = padj, n = 50) %>%
  pull(gene)

model_df <- as.data.frame(t(vsd_mat[selected_genes, ])) %>%
  rownames_to_column("sample") %>%
  left_join(metadata %>% rownames_to_column("sample"), by = "sample")

# Convert outcome to factor
model_df$condition <- factor(model_df$condition)

# Train/test split
set.seed(123)
train_index <- createDataPartition(model_df$condition, p = 0.8, list = FALSE)
train_data <- model_df[train_index, ]
test_data  <- model_df[-train_index, ]

# Logistic regression with regularization
x_train <- model.matrix(condition ~ . - sample, data = train_data)[, -1]
y_train <- train_data$condition

x_test <- model.matrix(condition ~ . - sample, data = test_data)[, -1]
y_test <- test_data$condition

cv_fit <- cv.glmnet(
  x = x_train,
  y = y_train,
  family = "binomial",
  alpha = 1
)

best_lambda <- cv_fit$lambda.min
best_lambda

pred_prob <- predict(cv_fit, newx = x_test, s = "lambda.min", type = "response")
pred_class <- ifelse(pred_prob > 0.5, levels(y_train)[2], levels(y_train)[1]) %>%
  factor(levels = levels(y_train))

confusionMatrix(pred_class, y_test)

# ROC curve
roc_obj <- roc(response = y_test, predictor = as.numeric(pred_prob))
auc(roc_obj)

plot(roc_obj, main = "ROC Curve for Biomarker-Based Classification")

# Random forest example
rf_model <- randomForest(
  condition ~ . - sample,
  data = train_data,
  importance = TRUE,
  ntree = 500
)

rf_pred <- predict(rf_model, newdata = test_data)
confusionMatrix(rf_pred, y_test)

# Variable importance
importance_df <- importance(rf_model) %>%
  as.data.frame() %>%
  rownames_to_column("feature") %>%
  arrange(desc(MeanDecreaseGini))

head(importance_df, 15)

ggplot(importance_df %>% slice_max(order_by = MeanDecreaseGini, n = 15),
       aes(x = reorder(feature, MeanDecreaseGini), y = MeanDecreaseGini)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Top Features by Random Forest Importance",
    x = "Feature",
    y = "Mean Decrease Gini"
  )

# Example: longitudinal biomarker measurements
biomarker_ts <- read.csv("biomarker_time_series.csv")

head(biomarker_ts)

# Suppose columns: date, patient_id, biomarker_value
biomarker_ts$date <- as.Date(biomarker_ts$date)

# Aggregate mean biomarker value by date
daily_signal <- biomarker_ts %>%
  group_by(date) %>%
  summarise(mean_value = mean(biomarker_value, na.rm = TRUE)) %>%
  arrange(date)

ggplot(daily_signal, aes(x = date, y = mean_value)) +
  geom_line() +
  labs(
    title = "Average Biomarker Signal Over Time",
    x = "Date",
    y = "Mean Biomarker Value"
  )

# Convert to time series object
signal_ts <- ts(daily_signal$mean_value, frequency = 7)

# Decomposition
signal_decomp <- stl(signal_ts, s.window = "periodic")
plot(signal_decomp)

# Automatic ARIMA model
fit_arima <- auto.arima(signal_ts)
summary(fit_arima)

# Forecast next 14 periods
signal_forecast <- forecast(fit_arima, h = 14)
plot(signal_forecast)

# Example metadata with repeated time points
metadata$timepoint <- factor(metadata$timepoint, levels = c("T0", "T1", "T2", "T3"))

dds_time <- DESeqDataSetFromMatrix(
  countData = counts_filtered,
  colData = metadata,
  design = ~ patient_id + timepoint
)

dds_time <- DESeq(dds_time)

# Compare T3 vs T0
res_time <- results(dds_time, contrast = c("timepoint", "T3", "T0"))
res_time_df <- as.data.frame(res_time) %>%
  rownames_to_column("gene") %>%
  filter(!is.na(padj)) %>%
  arrange(padj)

head(res_time_df)

# Plot trajectories for selected genes
trajectory_genes <- res_time_df %>%
  slice_min(order_by = padj, n = 6) %>%
  pull(gene)

traj_df <- vsd_mat[trajectory_genes, ] %>%
  as.data.frame() %>%
  rownames_to_column("gene") %>%
  pivot_longer(-gene, names_to = "sample", values_to = "expression") %>%
  left_join(metadata %>% rownames_to_column("sample"), by = "sample")

ggplot(traj_df, aes(x = timepoint, y = expression, group = patient_id, color = patient_id)) +
  geom_line(alpha = 0.7) +
  geom_point() +
  facet_wrap(~ gene, scales = "free_y") +
  labs(
    title = "Gene Expression Trajectories Across Time",
    x = "Time Point",
    y = "Variance-Stabilized Expression"
  )

# Convert gene symbols to ENTREZ IDs
gene_symbols <- sig_res$gene

gene_map <- bitr(
  gene_symbols,
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Hs.eg.db
)

head(gene_map)

# GO enrichment
ego <- enrichGO(
  gene = gene_map$ENTREZID,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.05,
  readable = TRUE
)

head(as.data.frame(ego))

# Dotplot
dotplot(ego, showCategory = 15)

# KEGG enrichment
ekegg <- enrichKEGG(
  gene = gene_map$ENTREZID,
  organism = "hsa",
  pvalueCutoff = 0.05
)

head(as.data.frame(ekegg))
barplot(ekegg, showCategory = 10)

# Save key results
write.csv(res_df, "deseq2_results_full.csv", row.names = FALSE)
write.csv(sig_res, "deseq2_significant_genes.csv", row.names = FALSE)
write.csv(importance_df, "model_feature_importance.csv", row.names = FALSE)

# Save transformed matrix
write.csv(as.data.frame(vsd_mat), "variance_stabilized_expression.csv")

# Session information
sessionInfo()

The post Digital Biology with R: Advanced Bioinformatics, Predictive Modeling, and Time Series Analysis for Modern Life Sciences appeared first on R Programming Books.

TLDR: To test whether I could tell the difference between decaf coffees I conducted a highly scientific test (subject to funding constraints).

One of my goals for 2025 was reducing my caffeine intake after having one too many sleepless nights. The problem was that all the decaffeinated coffee I’d tried was terrible.

Or was it?

After all, I once called myself an audiophile until a series of A/B tests suggested I couldn’t tell the difference between tracks encoded at different bitrates. So, it was entirely possible that I’d been brainwashed by big coffee to believe decaf coffee was inferior.

But, how exactly could we test this out?

The most obvious solution was to run a cross-country randomized double-blind experiment. This way, I wouldn’t automatically base my resentment on the caffeinated status of the coffee and focus on my subjective rating of the quality of each coffee.

Which is pretty much what I did:

Step 1: Sample selection

The first step was to select a wide enough sample of coffee beans to make the study as sciency as possible. Roping in my wife to help out, I purchased as many decaf varieties I could get our hands on.

Step 2: Sample blinding

After selecting a large representative sample of coffees (n=6), I packed a sample of each in its own container (pictured). To obscure each coffee’s identity I assigned them a number from 1 to 6. To further enhance the science I then had my wife assign new numbers so neither of us knew the origin of each sample.

Step 3: Testing

Before starting the test I cleaned and descaled the coffee machine. Beans from each container were freshly ground at room temperature and used to make six separate espressos. Shots were drawn on a quasi-random basis according to whatever my wife handed to me. We then took a sip of each coffee and ranked our preferences from one to six.

Step 4: Results

Although I’d have liked to pre-register my research, none of the top econometric journals I contacted expressed interest. However, my running assumption was that our preferences for a coffee were mainly psychological and had little to do with its caffeine content.

If this was true, I’d expect to see no relationship between our rankings. But, to my surprisethis didn’t appear to be the case. Instead, we both ranked the beans in a similar order:

Code Snippet:

#load libraries and import data
library(tidyverse)
dta_coffee_science<-read_csv("./Data/250216 blind coffee ratings.csv")

# Show linear association between samples by assigned label
#reverse axis so lower rankings are higher on the axis scale
plt_rankings_by_coffee_no <- ggplot(data = dta_coffee_science,
                                  aes(y = ranking_person_b, x = ranking_person_a)) +
 geom_text(aes(label = blind_label_round_2), size = 3.5) +
 scale_y_reverse(name = "Person A Ranking (1 = Best)") +
 scale_x_reverse(name = "Person B Ranking (1 = Best)") + 
 labs(title = "Coffee Sample Rankings: Person A vs Person B",
      subtitle = "Double-blind taste test results of brewed coffee samples",
      caption = "Note: Lower numbers indicate higher preference") +
 theme_classic()


plt_rankings_by_coffee_no

Of course, we’re doing some real science here, so to check let’s apply Kendall’s Tau of ranks and the Spearman correlation test for a null-hypothesis that there is no statistical association between our rankings.

Code Snippet:

#kendall
cor.test(data=dta_coffee_science, 
         ~ ranking_person_a + ranking_person_b, method = "kendall")

#spearman
cor.test(data=dta_coffee_science, 
         ~ ranking_person_a + ranking_person_b, method = "spearman")

With p-values from six to eight percent, this isn’t a ringing endorsement of the results, but having already written the blog I’m happy to adjust my definition of significant to conclude our preferences were similar to one another.

Of course, my willingness to play fast and lose with the stats also stems from knowing a key result: we both ranked the store-bought caffeinated beans highest.

Code Snippet:

# Show linear association between samples by assigned label and caffeination status 
plt_rankings_by_caffeine <-ggplot(data=dta_coffee_science,
                     aes(y=ranking_person_b, x=ranking_person_a,col=decaf))+
 geom_text(aes(label = blind_label_round_2), size = 3.5) +
  scale_y_reverse(limits = c(7.5, 0.5)) +   
  scale_x_reverse(limits = c(7.5, 0.5)) +  
    coord_cartesian(clip = "off") +  
 labs(title = "Coffee Sample Rankings: Person A vs Person B",
      subtitle = "Double-blind taste test results of brewed coffee samples",
      caption = "Note: Lower numbers indicate higher preference",
      x="Person A Ranking (1 = Best)", 
      y="Person B Ranking (1 = Best)") +
 theme_classic()+
 scale_color_manual(values = c("black", "blue"), name = "Decaf:") 

plt_rankings_by_caffeine

I also found it surprising that the beans from a specialized provider of decaf weren’t necessarily ranked higher, with only one of their beans ranked in the top three:

Code Snippet:

# Show linear association between samples by assigned label and caffeination status with original labels
plt_rankings_by_caffeine_named <- ggplot(data = dta_coffee_science,
                     aes(y = ranking_person_b, x = ranking_person_a, col = decaf)) +
  geom_text(aes(label = str_wrap(paste0(coffee_brand, ": ", coffee_name), width = 15)), 
            size = 3.5, lineheight = 0.8) +
  scale_y_reverse(limits = c(7.5, 0.5)) +   
  scale_x_reverse(limits = c(7.5, 0.5)) +   
  coord_cartesian(clip = "off") +            # stop clipping text at panel border
  labs(title = "Coffee Sample Rankings: Person A vs Person B",
       subtitle = "Double-blind taste test results of brewed coffee samples",
       caption = "Note: Lower numbers indicate higher preference",
       x = "Person A Ranking (1 = Best)", 
       y = "Person B Ranking (1 = Best)") +
  theme_classic() +
  theme(plot.margin = margin(10, 60, 10, 60)) +  
  scale_color_manual(values = c("black", "blue"), name = "Decaffeinated")

plt_rankings_by_caffeine_named

And while these coffee nerds might disagree, the results suggest we can tell the difference between coffees and both prefer the caffeinated alternative.

When I recounted the result to a food chemist they told me this probably has something to do with decaf coffee lacking the bitterness of caffeine.

When I recounted the results to my wife, she told me to never waste her time like this again. I probably will.

In the spirit of open science, you can download the dataset here.

The post Using science to find the best decaf appeared first on Giles.

Short practical advice on attribution:

1. Treat last-touch as a bias, not as a neutral baseline – If most advertisers optimize to last-touch conversions, budgets and auction pressure will drift toward channels that are easy to credit rather than channels that are truly incremental.

2. Advanced attribution is valuable because it changes spending behavior – The upside is not cleaner reporting. The upside is spending less on overpriced “credit-capturing” inventory and more on channels that create demand earlier in the journey.

3. Non-click measurement matters most in click-biased markets – If the market still rewards final clicks, then impression effects, sequence effects, and assisted conversions are systematically underpriced.

4. Do not expect more sophistication to automatically raise profit – Better measurement only helps if it changes bidding, budget allocation, and evaluation rules inside the company.

5. Add incrementality before adding complexity – A simple holdout, geo test, or lift study is often more useful than another dashboard built on the same click-based attribution logic.

Long Version

I just read Ron Berman’s paper Beyond the Last Touch: Attribution in Online Advertising, and I think the most useful implication for practitioners is slightly different from the headline result.

The paper compares no attribution, last-touch attribution, and Shapley-value attribution in a multi-publisher advertising market. Its central result is that attribution is not just a reporting layer. It changes bidding incentives and therefore changes market outcomes.

That matters because most advertisers still do not operate with advanced measurement. In many cases they rely on some mix of:

• last-touch or last-click attribution
• platform-reported conversions
• weak or infrequent incrementality testing
• budget allocation rules based on short-term ROAS

In that world, advanced attribution has a second source of value: it helps you avoid the distortions created by everyone else’s bad measurement.

What the paper actually shows

Berman models a market where one advertiser buys media across multiple publishers while other advertisers are more local. Consumers may see ads across different publishers, and one exposure can affect the value of another. That creates externalities across touchpoints.

The important result is that these externalities make optimization hard even when everyone is symmetric in information. This is not just a fraud problem or a platform-information problem. It is a structural feature of multi-publisher advertising.

The paper then shows two things that are highly relevant in practice:

1. Last-touch attribution often pushes advertisers to overbid.
2. A more balanced attribution rule such as the Shapley value usually performs better than last touch for the advertiser.

There is also a useful warning in the paper: more accurate attribution does not always increase advertiser profit. Once you account for market equilibrium and competitor response, the relationship between better measurement and better outcomes is not trivial.

That is the academic result. The practical extension is where things get interesting.

Why advanced attribution becomes more valuable when most advertisers still use last touch

If most of the market uses last-touch logic, then the market is not being priced on incremental contribution. It is being priced on who gets the final credit.

That sounds like a reporting issue, but it is really a budget allocation issue.

Channels that close demand tend to look stronger than they are. Branded search, retargeting, affiliate traffic, and other lower-funnel placements often appear to be doing all the work because they sit close to the conversion event. Channels that create demand earlier in the journey often look weaker because they rarely receive the final touch.

If enough advertisers optimize this way, a predictable pattern emerges:

• lower-funnel inventory becomes crowded and expensive
• channels that harvest existing intent absorb too much budget
• upper- and mid-funnel channels look worse than their true contribution
• advertisers confuse ease of measurement with causal impact

This is exactly where advanced attribution helps.

It helps not because it gives you a prettier customer journey chart, but because it gives you a better pricing model for media. If your attribution system captures assisted effects, non-click influence, or incremental lift more accurately than the market standard, you can avoid overpaying for touchpoints that are merely good at showing up last.

Why non-click measurement matters

The paper itself is about attribution rules, not a direct comparison of click-based versus non-click-based systems. But the implication is clear.

If the dominant market standard is click-biased last-touch measurement, then any method that can recover value from non-click exposures has a structural advantage.

That could include:

• impression-based attribution
• conversion lift studies
• geo experiments
• holdout tests
• media mix modeling
• sequence analysis using first-party exposure data

These methods are imperfect, and they answer slightly different questions. But they all do something last-touch usually cannot do well: they credit touchpoints that influence conversion without demanding a final click.

In a market full of click-based optimization, that means advanced measurement can uncover media that is genuinely incremental but systematically undervalued.

The real gain is not better reporting. It is better spending behavior.

I think this is the most important practical takeaway from the paper.

Advanced attribution should be judged by whether it improves decisions such as:

• Which channels deserve more budget?
• Which channels are just harvesting demand that already exists?
• Where are we overbidding because the market over-credits the final touch?
• Which publishers look weak only because our measurement ignores assisted effects?

If your attribution system changes those decisions, it has value.

If it only creates a more sophisticated dashboard while the company still optimizes to platform ROAS and last-click conversions, it has much less value.

That is also why many measurement projects disappoint. The model gets better, but the organization does not change the control system. Finance still trusts the old KPI. Paid media teams still chase the same targets. Creative teams still optimize for clicks because clicks are what get reported weekly.

The paper’s logic fits this well: attribution changes incentives. If the incentives do not change, the gains from attribution are limited.

A simple mental model for advertisers

If most competitors use last touch, assume three things:

1. Some inventory is overpriced because it captures credit, not because it creates lift.
2. Some inventory is underpriced because it influences conversion without closing it.
3. Your edge comes from estimating that gap better than the market.

That is why advanced attribution and incrementality measurement can be valuable even if they are noisy. You do not need a perfect model. You need a model that is less wrong than the one most of the market is using.

Where the paper is most useful, and where it is limited

One thing I like about the paper is that it is transparent about trade-offs. It does not claim that better attribution always improves everything. In fact, some attribution rules improve advertiser profit partly by shifting profit away from publishers, not by massively improving total market efficiency.

That is a useful correction to the usual industry narrative that “better measurement” always makes the whole system work better.

At the same time, the paper is stylized. It uses a game-theoretic model with a limited number of publishers and advertisers. It is not a direct operational guide to how to run MMM, lift testing, or conversion APIs in a modern stack.

So I would use it for the strategic lesson, not for a literal implementation recipe.

The strategic lesson is:

When the market overuses last-touch attribution, advanced measurement becomes more valuable because it helps advertisers resist the bidding and budgeting distortions created by everyone else’s simplistic measurement.

What I would do in practice

If I were working with an advertiser that still leans heavily on last-touch reporting, I would not start by building a highly complex attribution model.

I would start with a simpler sequence:

1. Keep last-touch reporting as an operational lens, not the source of truth.
2. Add one credible incrementality layer such as a holdout test, geo test, or lift study.
3. Separate channels that create demand from channels that mostly capture demand.
4. Use attribution outputs to inform budget ranges, not to automate every bid immediately.
5. Revisit evaluation metrics so teams are not punished for investing in channels that assist but rarely close.

That is probably the highest-return path for most advertisers because the biggest issue is usually not lack of modeling sophistication. It is overdependence on a biased attribution rule.

Bottom line

If most advertisers still use last touch and weak testing, advanced attribution is more valuable than it first appears.

Not because it makes reports look better, but because it can help you spend against incremental contribution while the rest of the market spends against credited contribution.

That gap is where the advantage is.

Paper: Ron Berman, Beyond the Last Touch: Attribution in Online Advertising.

Contributors: Michael Love, Justin Landis, Pierre-Paul Axisa

The generality of tidySummarizedExperiment makes it easy to interface with several tidyverse packages (e.g. dplyr, tidyr, ggplot2, purrr, plotly). This is possible thanks to its approach of converting SummarizedExperiment objects to tibbles, performing operations, and converting back to the original format. This conversion process introduces substantial overhead when working with large-scale datasets. Each operation requires multiple data transformations, with the conversion to tibble format creating memory copies of the entire dataset, followed by the reverse conversion back to SummarizedExperiment. For datasets containing hundreds of samples and tens of thousands of genes, these repeated conversions can consume memory and add significant computational overhead to even simple operations such as filtering or grouping.

With the new tidySummarizedExperiment release (v1.19.7), we have introduced new optimisations that address these performance limitations. This optimisation is powered by:

1. Check for the query domain (assay, colData, rowData), and execute specialised operation.
2. Use of plyxp for complex domain-specific queries.

plyxp is a tidyomics package developed by Justin Landis, and first released as part of Bioconductor 3.20 in October 2024. It uses data-masking functionality from the rlang package to perform efficient operations on SummarizedExperiment objects.

mutate.SummarizedExperiment <- function(.data, ...) {
    # Legacy implementation of mutate() for SummarizedExperiment:
    # - Validates requested edits against special/view-only columns
    # - Performs mutate() via tibble round-trip, then reconstructs the SE
    # Check that we are not modifying a key column
    cols <- enquos(...) |> names()
    
    # Deprecation of special column names:
    # capture all quoted args to detect deprecated special-column usage
    .cols <- enquos(..., .ignore_empty="all") %>% 
        map(~ quo_name(.x)) %>% unlist()
    if (is_sample_feature_deprecated_used(.data, .cols)) {
        # Record deprecated usage into metadata for backward compatibility
        .data <- ping_old_special_column_into_metadata(.data)
    }
    
    # Identify view-only/special columns (sample/feature keys, etc.)
    # Use a small slice to reduce overhead while probing structure
    special_columns <- get_special_columns(
        # Decrease the size of the dataset
        .data[1:min(100, nrow(.data)), 1:min(20, ncol(.data))]
    ) |> c(get_needed_columns(.data))
    
    # Are any requested targets among special/view-only columns?
    tst <-
        intersect(
            cols,
            special_columns
        ) |> 
        length() |>
        gt(0)

    if (tst) {
        columns <-
            special_columns |>
                paste(collapse=", ")
        stop(
            "tidySummarizedExperiment says:",
            " you are trying to rename a column that is view only",
            columns,
            "(it is not present in the colData).",
            " If you want to mutate a view-only column,",
            " make a copy and mutate that one."
        )
    }

    # If Ranges column not in query, prefer faster tibble conversion
    # Skip expanding GRanges columns when not referenced
    skip_GRanges <-
        get_GRanges_colnames() %in% 
        cols |>
        not()
    
    # Round-trip: SE -> tibble -> dplyr::mutate -> SE
    .data |>
        as_tibble(skip_GRanges=skip_GRanges) |>
        dplyr::mutate(...) |>
        update_SE_from_tibble(.data)
}

mutate.SummarizedExperiment <- function(.data, ...) {

       # Check if query is composed (multiple expressions)
    if (is_composed("mutate", ...)) return(decompose_tidy_operation("mutate", ...)(.data))

        # Check for scope and dispatch elegantly
        scope_report <- analyze_query_scope_mutate(.data, ...)
        scope <- scope_report$scope

        result <-
        if(scope == "coldata_only") modify_samples(.data, "mutate", ...)
        else if(scope == "rowdata_only") modify_features(.data, "mutate", ...)
        else if(scope == "assay_only") mutate_assay(.data, ...)
        else if(scope == "mixed") modify_se_plyxp(.data, "mutate", scope_report, ...)
        else mutate_via_tibble(.data, ...)

        # Record latest mutate scope into metadata for testing/introspection
        meta <- S4Vectors::metadata(result)
        if (is.null(meta)) meta <- list()
        meta$latest_mutate_scope_report <- scope_report
        S4Vectors::metadata(result) <- meta

        return(result)

}

suppressPackageStartupMessages({
  library(ggplot2)
  library(dplyr)
  library(SummarizedExperiment)
  library(rlang)
  library(devtools)
  library(airway)
  library(microbenchmark)
  library(reactable)
  library(patchwork)
})

load_branch_code <- function(worktree_dir) {
  if (!requireNamespace("devtools", quietly = TRUE)) stop("Please install devtools to run this vignette.")
  # Debug: print the directory we're looking for
  cat("Looking for worktree directory:", worktree_dir, "\n")
  cat("Directory exists:", dir.exists(worktree_dir), "\n")
  cat("Current working directory:", getwd(), "\n")
  # Check if directory exists
  if (!dir.exists(worktree_dir)) {
    stop(paste("Worktree directory does not exist:", worktree_dir))
  }
  suppressMessages(devtools::load_all(worktree_dir, quiet = TRUE))
}

create_airway_test_se <- function() {
  suppressPackageStartupMessages(library(airway))
  data(airway)
  se <- airway
  se[1:200, ]
}

benchmark_scenarios <- function() {
  list(
    coldata_simple_assignment = quo({ se |> mutate(new_dex = dex) }),
    coldata_arithmetic = quo({ se |> mutate(avgLength_plus_5 = avgLength + 5) }),
    coldata_concat = quo({ se |> mutate(sample_info = paste(cell, dex, SampleName, sep = "_")) }),
    coldata_grouped_mean = quo({ se |> group_by(dex) |> mutate(avgLength_group_mean = mean(avgLength)) |> ungroup() }),
    assay_simple_assignment = quo({ se |> mutate(counts_copy = counts) }),
    assay_plus_one = quo({ se |> mutate(counts_plus_1 = counts + 1) }),
    assay_log = quo({ se |> mutate(log_counts_manual = log2(counts + 1)) }),
    complex_conditional_coldata = quo({ se |> mutate(length_group = ifelse(avgLength > mean(avgLength), "longer", "shorter")) }),
    complex_nested = quo({ se |> mutate(complex_category = ifelse(dex == "trt" & avgLength > mean(avgLength), "treated_long", ifelse(dex == "untrt", "untreated", "other"))) }),
    mixed_assay_coldata = quo({ se |> mutate(new_counts = counts * avgLength) }),
    multiple_simple_assay = quo({ se |> mutate(normalized_counts = counts / 1000, sqrt_counts = sqrt(counts)) }),
    chained_mutates = quo({ se |> mutate(tmp = avgLength * 2) |> mutate(flag = ifelse(tmp > mean(tmp), 1, 0)) }),

    # Filter benchmarks (scoped and non-rectangular)
    filter_coldata_simple = quo({ se |> filter(dex == "trt") }),
    filter_coldata_numeric = quo({ se |> filter(avgLength > median(avgLength)) }),
    filter_assay_nonrect = quo({ se |> filter(counts > 0) }),

    # Select benchmarks (covering colData-only, rowData-only, assays-only, mixed)
    select_coldata_simple = quo({ se |> select(.sample, dex) }),
    select_rowdata_simple = quo({ se |> select(.feature) }),
    select_assay_only = quo({ se |> select(counts) }),
    select_mixed_keys_counts = quo({ se |> select(.sample, .feature, counts) }),
    select_coldata_wide = quo({ se |> select(.sample, dex, avgLength, SampleName) }),

    # Distinct benchmarks (covering colData-only, rowData-only, assays-only, mixed)
    distinct_coldata_simple = quo({ se |> distinct(dex) }),
    distinct_coldata_multiple = quo({ se |> distinct(dex, avgLength) }),
    distinct_rowdata_simple = quo({ se |> distinct(.feature) }),
    distinct_assay_only = quo({ se |> distinct(counts) }),
    distinct_mixed_keys_counts = quo({ se |> distinct(.sample, .feature, counts) }),
    distinct_coldata_wide = quo({ se |> distinct(.sample, dex, avgLength, SampleName) }),
    distinct_with_keep_all = quo({ se |> distinct(dex, .keep_all = TRUE) }),
    distinct_complex_expression = quo({ se |> distinct(dex, avgLength) })
  )
}

run_one <- function(expr_quo, reps = 5L) {
  se_base <- create_airway_test_se()
  mb <- microbenchmark::microbenchmark(
    eval_tidy(expr_quo),
    times = reps,
    setup = { se <- se_base },          # reuse the same input, avoid recreating inside the timed expr
    control = list(warmup = 2L)
  )
  # microbenchmark returns nanoseconds; convert to milliseconds
  as.numeric(mb$time) / 1e6
}

run_all_scenarios <- function(branch_label, reps = 7L) {
  scenarios <- benchmark_scenarios()
  out <- list()
  for (nm in names(scenarios)) {
    tms <- run_one(scenarios[[nm]], reps = reps)
    out[[length(out) + 1L]] <- data.frame(
      branch = branch_label,
      scenario = nm,
      replicate = seq_along(tms),
      elapsed_ms = tms,
      stringsAsFactors = FALSE
    )
  }
  bind_rows(out)
}

# Parallel version: run each scenario on a separate worker
run_all_scenarios_parallel <- function(branch_label, reps = 20L, workers = 1L, initializer = NULL) {
  scenarios <- benchmark_scenarios()
  nms <- names(scenarios)
  old_plan <- future::plan()
  on.exit(future::plan(old_plan), add = TRUE)
  future::plan(future::multisession, workers = workers)
  res <- future.apply::future_lapply(nms, function(nm) {
    if (!is.null(initializer)) initializer()
    tms <- run_one(scenarios[[nm]], reps = reps)
    data.frame(
      branch = branch_label,
      scenario = nm,
      replicate = seq_along(tms),
      elapsed_ms = tms,
      stringsAsFactors = FALSE
    )
  }, future.seed = TRUE)
  dplyr::bind_rows(res)
}

# Worktree directories (already exist in the post directory)
wt_before <- ".__bench_before__"
wt_after <- ".__bench_after__"

# Verify worktrees exist
if (!dir.exists(wt_before)) {
  stop("Worktree directory does not exist: ", wt_before)
}
if (!dir.exists(wt_after)) {
  stop("Worktree directory does not exist: ", wt_after)
}

# Before optimisation (commit 87445757)
load_branch_code(wt_before)

Looking for worktree directory: .__bench_before__ 
Directory exists: TRUE 
Current working directory: /Users/a1234450/Documents/GitHub/tidyomicsBlog/posts/2025-10-25-tidySummarizedExperiment-optimization 

res_before <- run_all_scenarios(branch_label = "before_optimization", reps = 10L)

# After optimisation (commit 9f7c26e)
load_branch_code(wt_after)

Looking for worktree directory: .__bench_after__ 
Directory exists: TRUE 
Current working directory: /Users/a1234450/Documents/GitHub/tidyomicsBlog/posts/2025-10-25-tidySummarizedExperiment-optimization 

res_after <- run_all_scenarios(branch_label = "after_optimization", reps = 10L)

results <- dplyr::bind_rows(res_before, res_after) |>
  dplyr::mutate(operation = dplyr::case_when(
    grepl("^filter", scenario) ~ "filter",
    grepl("^select", scenario) ~ "select",
    grepl("^distinct", scenario) ~ "distinct",
    TRUE ~ "mutate"
  ))

summary_table <- results |>
  group_by(branch, scenario) |>
  summarise(median_ms = median(elapsed_ms), .groups = "drop") |>
  tidyr::pivot_wider(names_from = branch, values_from = median_ms) |> 
  dplyr::mutate(speedup = round(before_optimization / after_optimization, 2))

dodge_w <- 0.7

p_box <- ggplot(results, aes(x = scenario, y = elapsed_ms, fill = branch)) +
  geom_boxplot(position = position_dodge(width = dodge_w), width = 0.7, outlier.shape = NA) +

  # Add jittered points aligned with the dodged boxplots
  geom_point(
    position = position_jitterdodge(jitter.width = 0.1, jitter.height = 0, dodge.width = dodge_w), 
    alpha = 0.6, 
    size = 0.5
  ) +
  scale_y_log10() + 
  coord_flip() +
  facet_grid(operation ~ ., scales = "free_y", space = "free_y") +
  annotation_logticks(sides = "b") +
  labs(title = "Performance comparison: Before vs After optimization",
       x = "Scenario",
       y = "Elapsed (ms)") +
  theme_bw() +
  
  # Angle x labels  
  theme(legend.position = "top", axis.text.x = element_text(angle = 45, hjust = 1))

# Speedup summary panel (median before/after ratio)
speedup_plot_data <- summary_table |>
  dplyr::mutate(operation = dplyr::case_when(
    grepl("^filter", scenario) ~ "filter",
    grepl("^select", scenario) ~ "select",
    grepl("^distinct", scenario) ~ "distinct",
    TRUE ~ "mutate"
  ))

p_speedup <- ggplot(
  speedup_plot_data,
  aes(x = speedup, y = reorder(scenario, speedup))
) +
  geom_col(width = 0.7, fill = "grey70", color = "grey40") +
  facet_grid(operation ~ ., scales = "free_y", space = "free_y") +
  labs(
    title = "Median speedup by scenario",
    x = "Speedup (before/after, x)",
    y = NULL
  ) +
  theme_bw() +
  theme(legend.position = "none")

combined_plot <- p_box + p_speedup + patchwork::plot_layout(widths = c(2.3, 1))
combined_plot

# Save the combined figure
ggsave("benchmark_plot.png", plot = combined_plot, width = 14, height = 8)

R version 4.5.3 (2026-03-11)
Platform: x86_64-apple-darwin20
Running under: macOS Sonoma 14.6.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Adelaide
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] tidySummarizedExperiment_1.19.7 tidyr_1.3.2                    
 [3] testthat_3.3.2                  ttservice_0.5.3                
 [5] patchwork_1.3.2                 reactable_0.4.5                
 [7] rlang_1.1.7                     microbenchmark_1.5.0           
 [9] airway_1.30.0                   SummarizedExperiment_1.40.0    
[11] Biobase_2.70.0                  GenomicRanges_1.62.1           
[13] Seqinfo_1.0.0                   IRanges_2.44.0                 
[15] S4Vectors_0.48.0                BiocGenerics_0.56.0            
[17] generics_0.1.4                  MatrixGenerics_1.22.0          
[19] matrixStats_1.5.0               dplyr_1.2.0                    
[21] ggplot2_4.0.2                   devtools_2.4.6                 
[23] usethis_3.2.1                  

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1    viridisLite_0.4.3   farver_2.1.2       
 [4] S7_0.2.1            fastmap_1.2.0       lazyeval_0.2.2     
 [7] digest_0.6.39       plyxp_1.4.3         lifecycle_1.0.5    
[10] ellipsis_0.3.2      magrittr_2.0.4      compiler_4.5.3     
[13] tools_4.5.3         yaml_2.3.12         data.table_1.18.2.1
[16] knitr_1.51          S4Arrays_1.10.1     labeling_0.4.3     
[19] htmlwidgets_1.6.4   pkgbuild_1.4.8      DelayedArray_0.36.0
[22] RColorBrewer_1.1-3  pkgload_1.5.0       abind_1.4-8        
[25] withr_3.0.2         purrr_1.2.1         desc_1.4.3         
[28] grid_4.5.3          fansi_1.0.7         scales_1.4.0       
[31] cli_3.6.5           rmarkdown_2.30      ragg_1.5.1         
[34] remotes_2.5.0       otel_0.2.0          rstudioapi_0.18.0  
[37] httr_1.4.8          sessioninfo_1.2.3   cachem_1.1.0       
[40] stringr_1.6.0       XVector_0.50.0      vctrs_0.7.1        
[43] Matrix_1.7-4        jsonlite_2.0.0      systemfonts_1.3.2  
[46] crosstalk_1.2.2     plotly_4.12.0       glue_1.8.0         
[49] reactR_0.6.1        stringi_1.8.7       gtable_0.3.6       
[52] tibble_3.3.1        pillar_1.11.1       htmltools_0.5.9    
[55] brio_1.1.5          R6_2.6.1            textshaping_1.0.5  
[58] rprojroot_2.1.1     evaluate_1.0.5      lattice_0.22-9     
[61] memoise_2.0.1       SparseArray_1.10.9  xfun_0.56          
[64] fs_1.6.7            pkgconfig_2.0.3    

© 2025 tidyomics. Content is published under Creative Commons CC-BY-4.0 License for the text and BSD 3-Clause License for any code. | R-Bloggers

This post is the last in a series of seven on population and people movement issues in the Pacific. The first six posts featured code re-generating the charts I used in a keynote speech before the November 2025 meeting of the Pacific Heads of Planning and Statistics in Wellington, New Zealand. Today’s post is simply a narrative drawing on all those charts. There’s no R code today; links to the previous posts are at the bottom.

Population growth is varied, but the larger Pacific island countries are growing pretty fast

The first point to make is that the Pacific is very varied in terms of its population dynamics. We can see this a bit in the most straightforward and intuitive chart of population historical growth and future projections:

The data this comes from is the United Nations’ 2024 population projections, which are currently used as-is in the Pacific Data Hub. Some of these figures are known to be wrong, either because there have been censuses that reported too late to be taken into account (e.g. Federated States of Micronesia) or for other reasons (e.g. Tokelau). But they’re a good starting point.

The vertical scale on that chart is “free” meaning each country-facet is on a different scale, so it’s not easy to tell visually which countries are larger or smaller. To help with this, they are organised from smallest (Niue) to largest (Papua New Guinea). Pitcairn Islands, although a member of the Pacific Community, is not included because of its tiny size (around 50 people) even by Pacific standards.

A bit of familiarity with the region though, helped by the sequencing on the chart, identifies that the larger countries—those in the bottom row from Samoa to Papua New Guinea, up to say Kiribati and Guam in the row above—are growing faster as well as being larger. Well that’s hardly surprising is it; faster growing countries will of course become larger, so we’d expect these things to be related? Except that the causality is not all in one direction like that. At its extreme, there’s more to the fact that Niue is only one or two thousand people and Papua New Guinea is 10,000 times its size (about 11 million) than that Papua New Guinea has grown faster recently. Niue’s small physical land mass, and other characteristics, are one of the reasons it doesn’t grow as fast.

The next chart is an attempt to illustrate this further. First, consider a version where we just put size on the horizontal axis and growth on the vertical. Both axes are logarithmically transformed, or else the chart would be all blank space with Papua New Guinea out on the right and all other countries a cluster of dots on top of eachother on the left.

I quite like this chart, for audiences with enough numeracy to cope with a scatter plot and log axes. The pink shading makes clear which countries are shrinking, a point of definite interest. The colour-coding of points by sub-region is very useful, making clear how the Melanesian countries cluster together in the top right “large and growing fast” quadrant; most Polynesian countries (not all) are shrinking; and Micronesia has a real mix. The next version of this chart adds a particularly interesting element to this (in a slide show, you can move from one chart to the next and it looks like these circles appear by magic):

Those circles represent the countries or territories that have some kind of easy migration access to a larger, richer country. This includes (conveniently, mostly in sets of three):

• three French territories (French Polynesia, Wallis and Futuna, New Caledonia);
• three Realm of New Zealand territories (Cook Islands, Tokelau, Niue);
• Pitcairn Island whose residents (mostly) have right of abode in the UK;
• three USA territories (Guam, Commonwealth of the Northern Mariana Islands, American Samoa)
• three independent countries with a Compact of Free Association with the USA (Palau, Marshall Islands, Federated States of Micronesia)

Now we have an interesting feature. All of the countries with negative growth apart from Tonga and Tuvalu are highlighted this way. Perhaps in future we would include Tuvalu in the list above because it does now have special arrangements with Australia and New Zealand that allow a certain amount of people movement, most importantly the Falepili Mobility Pathway scheme that allows 280 permanent residence visas per year with Australia. But this is too recent to show up in the chart above.

Samoa and Tonga do not have automatic right of residence to New Zealand but they have strong cultural and historical ties, large communities already living there (more on this later), and the visa obstacles are largely surmountable.

In my view the distinction between the countries that have easy mobility to a larger richer countries, and those that don’t, is the most important single marker to use when considering population issues in the Pacific.

Net migration is a critical factor in different population growth rates

OK, let’s look at net migration. The United Nations has to estimate this as part of its projection process so we can see it for all the countries in their dataset. In the below, I’ve organised countries from those least impacted by migration (top left) to most (bottom right).

There’s a few remarkable things here. One is that the chart is mostly red—migration is nearly all outwards. Another is that it leaps around a bit. We can suspect data problems for some of this eg in the case of Tokelau’s recent years, I just think we have an error (it’s being looked at).

As a way of drilling into the impact of migration on population, for just a subset of countries because it’s getting complicated, is comparing the natural rate of increase (ie births minus deaths) with the net migration—which between them, add up to total change in population for each country. Here’s a chart that does that, for just six countries:

In each of these charts, we can start with the green dashed line, which is where natural increase would be with no migration. Then we add (or subtract, in most cases) the red dotted line which is net migration. The sum of these two is the blue solid line, total population change. The sort of things we see here are:

• For Kiribati and Papua New Guinea, the solid blue line is fairly close to the dashed green line, indicating that natural increase is what is driving population change
• For Samoa and Marshall Islands, the green dashed natural increase line is very positive and the red dotted migration line is very negative. These two drivers mostly cancel eachother out, but when all is added up Marshall Islands has rapid recent population decline and Samoa still has some remnant slow population growth
• For Niue and Marianas, the story is more complicated but in recent years has stabilised at “not much change”

What we’re seeing here is that migration—or the lack of it in significant numbers, for Kiribati and Papua New Guinea—is what is driving the population story.

Migration (short and long term) changes the shape of the origin country’s demographics

What does this mean for the structure of who is left? To illustrate this, I like to compare Kiribati and Marshall Islands. Both are entirely or mostly coral atolls; they are only around an hour’s flight from eachother; and they have few natural resources other than their people, the ocean and its fish, and location.

Location is a critical asset or curse for the Marshall Islands. Kwajalein Atoll in the Marshalls was a major Japanese base in World War II and the site of a bloody battle in 1944; now it is a key US base forming a bridging zone and operational depth between the so-called second (Guam, Palau, Saipan, etc.) and third (Aleutians, Hawaii, Samoa, New Zealand) island chains in preparation for the next Pacific war, against whichever eastwards-facing Asian land power that might be. Bikini Atoll in the Marshalls was the site of US nuclear weapon testing during the cold war. For our purposes, all this matters because Marshall Islands has a Compact of Free Association with the USA which provides large amounts of funding plus free people movement to the USA. Even with recent crackdowns in the USA this persists, although dealing with Marshallese who have been forcibly returned from the USA due to low level criminal behaviour is becoming a policy challenge.

Kiribati has no such arrangement with its own large regional partner, Australia (former-British Kiribati uses the Australian dollar and drives on the left, just as Marshall Islands uses the US dollar and drives on the right hand side of the road).

To repeat some text from the blog post where I introduced this chart: Kiribati today has about four times the population of Marshall Islands but in 1980 was only about double. The significant thing here is the wasp waist of the Marshall Islands pyramid in 2025—while it had a similar shape to Kiribati in 1980. People at peak working and reproductive age are literally absent from today’s Marshall Islands—in this case, primarily in the USA.

The result of this is that Marshall Islands not only benefits from its individuals having more freedom of movement and opportunity, and sending back remittances from relatively high paying lives in the USA; but also having a pressure valve for what would otherwise be a rapidly (too fast?) growing population. To put it bluntly, Kiribati has a problem of too many people (particularly on crowded southern Tarawa); Marshall Islands, if it has a population problem, is one of too few. The contrast of crowded, relatively poor Tarawa and less-crowded, relatively well-off Majuro is an obvious and stark one to anyone travelling to them both in quick succession.

Pacific people are overseas in very considerable numbers

OK, so people have been moving from the Pacific islands to elsewhere for decades or longer. Proportionately speaking, does this matter? Have large numbers of Pacific islanders cumulatively ended up elsewhere? The following chart answers this with a resounding “yes”, for at least nine countries:

For seven countries, there are more people ethnically associated with that country living in the USA, New Zealand and Australia than in the origin country. For the countries in the bottom row – the three New Zealand Realm countries plus Pitcairn – there are many multiples more people living overseas. Around 40,000 Niueans live overseas (mostly in New Zealand) and less than 2,000 in Niue. For Tonga, Samoa and Marshall Islands the situation is not as extreme but still very substantial.

Illustrating this further, consider this chart of the world’s largest Pacific Islander cities. Thanks to comments on LinkedIn I have improved this from previous versions—we now have a better estimate of people with Pacific Islander as one ‘race’ of several in Hawaii, and a more comparable definition of Greater Wellington. But what we see is broadly the same message as before—of the top ten Pacific Islander cities in the world, two are in Australia (Sydney and Brisbane) and one each in New Zealand (Auckland) and the USA (Honolulu).

All this has substantial implications for cultural and national identity, for economics and for politics. What political rights and expectations are there for these overseas people? For the home country, is there a minimum size that is viable? Or even a minimum proportion of your people? What happens legally, politically and diplomatically (for example, to votes in the United Nations General Assembly) if—as is clearly possible under climate change—a low-lying coral atoll country like Tokelau, Marshall Islands or Tuvalu loses all of its land to sea-level rise and all of its people are living overseas?

Remittances are a critical, even dominant, part of many (but not all) Pacific island economies

Without going to these extreme scenarios we have enormous economic implications in the here and now. One of the strongest indicators of this is the level of remittances. Remittances are payments from family or other contacts overseas, typically in a higher income country. The source of remittances can be people on relatively short trips overseas—in the Pacific, examples include people in the Pacific Australia Labour Mobility scheme or the New Zealand Recognised Seasonal Employer scheme—or from long term migrants who have made the other country their indefinite home.

We see from this, final, chart that some Pacific island countries have extraordinarily high levels of remittances compared to averages of comparable countries, including other small states:

Interestingly, the three highest countries on this measure are not on my list of countries with special access to permanent residency in a large rich country. But as previously mentioned, Samoa and Tonga have particularly strong ties to New Zealand that act partway to such special access. Vanuatu does not, but it is the beneficiary of short term labour schemes. It’s also possible that remittances are under-reported in some countries.

For countries like Tonga, it seems likely or at least possible that remittances are coming from long term migrants in New Zealand. After a generation or three, will they stop sending money back to Tonga? Possibly, but as we have seen there is a continual ongoing refreshment in the form of new migrants and it is likely the remittances will continue for the forseeable future.

One thing is clear—labour mobility is an enduring feature of the Pacific region that meets multiple groups’ needs and leaves a decisive mark on both sending and receiving economies. Remittances are just one most direct part of the economic impact; others include investment ties and human capability development. But I won’t go into a literature review on this area.

There are important policy implications

Why do we care? There are important implications of all this. The context of my original talk was a meeting bringing together heads of National Statistical Offices with heads of national planning, and I wanted to highlight population issues as one area where the links between official statistics and national planning and policy are (or should be) particularly strong. Some of the key economic and planning issues in this regard as I see them include:

• Population projections (including at sub-national level) and the obvious implications for infrastructure and related planning (if and where to build roads, medical centres, schools, etc.)
• Impacts of migration (short and long term) on the working-age population back at home and what planning or policy levers are needed to deal with this, not just lament it
• Social impacts of many working-age people, perhaps disproportionately of one gender, being overseas
• Implications for taxation policy of a large proportion of national income coming in remittances
• Impact of migration (in either direction) on the net fiscal position (is it working age people coming to our country, paying more tax than they extract in benefits— or vice versa)
• Specifically, where do retirees end up, and who pays for them
• Impact of returnees from overseas work experiences—including raised skills and experiences, higher taxation payments, and general impact on national capabilities and capacity
• Impacts on equality back home from different types of people in the diaspora
• What labour market and cultural context should education be preparing young people for?

This aims to be an indicative, rather than comprehensive, list. After all, my aim here is to highlight some issues, not venture into policy advice that I have no mandate for and which would require much more systematic evaluation of options and root causes of the problems.

Other posts in this series

The seven blog posts in total in this series are set out below. The first six contain R code and data sources for each chart:

• Visual summaries of population size and growth
• Net migration
• Population pyramids
• World cities with the most Pacific Islanders
• Pacific diaspora
• Remittances
• Tying it all together (this post today)